{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# ML_in_Finance-Bayesian_vs_Frequentist_Estimation\n",
    "# Version: 1.0 (28.4.2020)\n",
    "# License: MIT\n",
    "# Email: paul@thalesians.co.uk and matthew.dixon@iit.edu\n",
    "# Notes: tested on Mac OS X running Python 3.6.9 with the following packages:\n",
    "# scipy=1.4,1, numpy=1.18.1, matplotlib=3.1.3, pandas=1.0.3\n",
    "# Citation: Please cite the following reference if this notebook is used for research purposes:\n",
    "# Dixon M.F., Halperin I. and P. Bilokon, Machine Learning in Finance: From Theory to Practice, Springer Graduate Textbook Series, 2020."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview\n",
    "The purpose of this notebook is to illustrate the difference between Frequentist and Bayesian estimation under Bernoulli trials. Please refer to Sections 3 and 4 of Chapter 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequentist versus Bayesian estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A random experiment: a single coin toss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider an experiment consisting in a single coin flip. We set the random variable $X$ to 0 if tails come up and 1 if heads come up. Then the probability density of $X$ is given by\n",
    "$$p(x \\, | \\, \\theta) = \\theta^x (1 - \\theta)^{1 - x},$$\n",
    "where $0 \\leq \\theta \\leq 1$ is the probability of heads showing up."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We view $p$ as a function of $x$, but parameterised by the given parameter $\\theta$, hence the notation, $p(x \\,|\\, \\theta)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will recognise $X$ as a **Bernoulli random variable**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A random experiment: multiple independent coin tosses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generalising somewhat, suppose that we perform $n$ such independent experiments (tosses) on the same coin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now,\n",
    "$$\\mathbf{X} = \\begin{pmatrix} X_1 \\\\ X_2 \\\\ \\vdots \\\\ X_n \\end{pmatrix} \\in \\mathbb{R}^n,$$\n",
    "where, for $1 \\leq i \\leq n$, $X_i$ is the result of the $i$th toss."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the probability density of $\\mathbf{X}$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the coin tosses are independent, the probability density of $\\mathbf{X}$, i.e. the joint probability density of $X_1, X_2, \\ldots, X_n$, is given by the product rule\n",
    "$$p(\\mathbf{x} \\,|\\, \\theta) = p(x_1, x_2, \\ldots, x_n \\,|\\, \\theta) = \\prod_{i=1}^n \\theta^{x_i} (1 - \\theta)^{1 - x_i} = \\theta^{\\sum x_i} (1 - \\theta)^{n - \\sum x_i}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose that we have actually observed ten coin tosses. On seven of them heads have come up, on the remaining three, tails. For example, we could have something like\n",
    "$$x_1 = 1, x_2 = 0, x_3 = 1, x_4 = 1, x_5 = 1, x_6 = 0, x_7 = 1, x_8 = 0, x_9 = 1, x_{10} = 1,$$\n",
    "although we shall see that the actual order in which heads and tails came up is unimportant. What really matters is that $n = 10$ and $\\sum_{i = 1}^n x_i = 7$.\n",
    "\n",
    "How can we use this information to **infer** or **estimate** the parameter $\\theta$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Frequentist estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Frequentist estimation: maximum likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both the frequentists and Bayesians regard the density $p(\\mathbf{x} \\,|\\, \\theta)$ as **likelihood**. Bayesians stick with this notation, whereas frequentists reinterpret $p(\\mathbf{x} \\,|\\, \\theta)$, which is a function of $\\mathbf{x}$ (given the parameters $\\mathbf{\\theta}$: in our case, there is a single parameter, so $\\theta$ is univariate, but this doesn't have to be the case) as a function of $\\theta$ (given the specific sample $\\mathbf{x}$), and write\n",
    "$$\\mathcal{L}(\\theta) := \\mathcal{L}(\\theta \\,|\\, \\mathbf{x}) := p(\\mathbf{x} \\,|\\, \\theta).$$\n",
    "Notice that we have merely reinterpreted this probability density, whereas its functional form remains the same, in our case:\n",
    "$$\\mathcal{L}(\\theta) = \\theta^{\\sum x_i} (1 - \\theta)^{n - \\sum x_i}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likelihood is one of the key ideas of the frequentist school. It was introduced by one of its founding fathers, Sir Ronald Aylmer Fisher:\n",
    "\n",
    " **\"*What has now appeared is that the mathematical concept of probability is ... inadequate to express our mental confidence or [lack of confidence] in making ... inferences, and that the mathematical quantity which usually appears to be appropriate for measuring our order of preference among different possible populations does not in fact obey the laws of probability. To distinguish it from probability, I have used the term \"likelihood\" to designate this quantity...*\"** &mdash; R.A. Fisher, *Statistical Methods for Research Workers*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is generally more convenient to work with the log of likelihood &mdash; the **log-likelihood**. Since $\\ln$ is a monotonically increasing function of its argument, the same values of $\\theta$ maximise the log-likelihood as the ones that maximise the likelihood.\n",
    "$$\\ln \\mathcal{L}(\\theta) = \\ln \\left\\{ \\theta^{\\sum x_i} (1 - \\theta)^{n - \\sum x_i} \\right\\} = \\left(\\sum x_i \\right) \\ln \\theta + \\left(n - \\sum x_i\\right) \\ln(1 - \\theta).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we find the value of $\\theta$ that maximises this expression? As in school calculus, we differentiate with respect to theta and solve for $\\theta$ that sets the (partial) derivative to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial}{\\partial \\theta} \\ln \\mathcal{L}(\\theta) = \\frac{\\sum x_i}{\\theta} + \\frac{n - \\sum x_i}{\\theta - 1}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equating this to zero and solving for $\\theta$, we obtain the **maximum likelihood estimate** for $\\theta$:\n",
    "$$\\hat{\\theta}_{\\text{ML}} = \\frac{\\sum x_i}{n}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To confirm that this value does indeed *maximise* the log-likelihood, we take the second derivative with respect to $\\theta$,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial^2}{\\partial \\theta^2} \\ln \\mathcal{L}(\\theta) = -\\frac{\\sum x_i}{\\theta^2} - \\frac{n - \\sum x_i}{(\\theta - 1)^2} < 0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this quantity is strictly negative for all $0 \\leq \\theta \\leq 1$, it is negative at $\\hat{\\theta}_{\\text{ML}}$, and we do indeed have a maximum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $\\hat{\\theta}_{\\text{ML}}$ depends only on the sum of the $x_i$s, we can answer our question: if in a sequence of 10 coin tosses exactly seven heads come up, then\n",
    "$$\\hat{\\theta}_{\\text{ML}} = \\frac{\\sum x_i}{n} = \\frac{7}{10} = 0.7.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we end up with a *single* value (a single \"point\") as our estimate, 0.7, in this sense we are doing **point estimation**. When we apply a Bayesian approach to the same problem, we shall see that the Bayesian estimate is a probability distribution, rather than a single point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have done quite a lot of mathematical work, but the answer is intuitively obvious. If we toss a coin ten times, and out of those ten times it lands with heads up, it is natural to estimate the probability of getting heads as 0.7. It's encouraging that the result of our maths agrees with our intuition and common sense."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assessing the quality of our estimator: bias and variance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we obtained our maximum likelihood estimate, we plugged in a specific number for $\\sum x_i$, 7, in this sense the estimator is an ordinary function.\n",
    "\n",
    "However, we could also view it as a function of the *random* sample,\n",
    "$$\\hat{\\theta}_{\\text{ML}} = \\frac{\\sum X_i}{n} = \\frac{7}{10} = 0.7,$$\n",
    "each $X_i$ being a random variable. A function of random variables is itself a random variable, so we can compute its expectation and variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In particular, an expectation of the **error**\n",
    "$$\\mathbf{e} = \\hat{\\mathbf{\\theta}} - \\mathbf{\\theta}$$\n",
    "is known as **bias**,\n",
    "$$\\text{bias}(\\hat{\\mathbf{\\theta}}, \\mathbf{\\theta}) = \\mathbb{E}(\\mathbf{e}) = \\mathbf{E}\\left[\\hat{\\mathbf{\\theta}}, \\mathbf{\\theta}\\right] = \\mathbf{E}\\left[\\hat{\\mathbf{\\theta}}\\right] - \\mathbf{E}\\left[\\mathbf{\\theta}\\right].$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As frequentists, we view the true value of $\\theta$ as a single, deterministic, fixed point, so we take it outside of the expectation:\n",
    "$$\\text{bias}(\\hat{\\mathbf{\\theta}}, \\mathbf{\\theta}) = \\mathbf{E}\\left[\\hat{\\mathbf{\\theta}}\\right] - \\mathbf{\\theta}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case it is\n",
    "$$\\mathbb{E}[\\hat{\\theta}_{\\text{ML}} - \\theta] = \\mathbb{E}[\\hat{\\theta}_{\\text{ML}}] - \\theta = \\mathbb{E}\\left[\\frac{\\sum X_i}{n}\\right] - \\theta = \\frac{1}{n} \\sum \\mathbb{E}[X_i] - \\theta = \\frac{1}{n} \\cdot n(\\theta \\cdot 1 + (1 - \\theta) \\cdot 0) - \\theta = 0,$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the bias is zero, so this particular maximum likelihood estimator is **unbiased** (otherwise it would be **biased**)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about the variance of this estimator?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\text{Var}[\\hat{\\theta}_{\\text{ML}}] = \\text{Var}\\left[\\frac{\\sum X_i}{n}\\right] \\overset{\\text{independence}}{=} \\frac{1}{n^2} \\sum \\text{Var}[X_i] = \\frac{1}{n^2} \\cdot n \\cdot \\theta (1 - \\theta) = \\frac{1}{n}\\theta(1 - \\theta),$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and we see that the variance of the estimator depends on the *true* value of $\\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, let us sample $n\\hat{\\theta}_{ML} \\sim Ber(\\theta, n)$ under the assumption of a biased coin. If we assume the number of Bernoulli trials is $n=10$ and then increase the sample size, $M$, of $\\hat{\\theta}_{ML}$ we can see that the sample mean of $\\hat{\\theta}_{ML}$ convergences to $\\theta$ (i.e. the sample bias of the estimator is zero) and the sample variance of $\\hat{\\theta}_{ML}$ converges to the true variance of the estimator $\\frac{1}{n}\\theta(1-\\theta)=0.021$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       M  sample bias  sample variance  true variance\n",
      "      10     1.00e-02           0.0166         0.0210\n",
      "     100     9.00e-03           0.0224         0.0210\n",
      "    1000    -3.10e-03           0.0219         0.0210\n",
      "   10000    -1.44e-03           0.0211         0.0210\n",
      "  100000     1.60e-04           0.0209         0.0210\n",
      " 1000000     2.09e-04           0.0210         0.0210\n",
      "10000000     5.54e-05           0.0210         0.0210\n"
     ]
    }
   ],
   "source": [
    "theta = 0.7\n",
    "n = 10 # the length of the Bernoulli sequence\n",
    "print(\"       M  sample bias  sample variance  true variance\")\n",
    "\n",
    "# M is the number of samples of n*theta ~ Ber(theta, n)\n",
    "for M in [10, 100, 1000, 10000, 100000, 1000000, 10000000]: \n",
    "    sample_theta = np.random.binomial(n, theta, size=M) / n\n",
    "    sample_bias_theta = np.mean(sample_theta) - theta\n",
    "    sample_var_theta = np.sum((sample_theta - np.mean(sample_theta))**2) / (M-1)\n",
    "    true_var_theta = theta * (1 - theta) / n\n",
    "    print(\"{:8} {: 12.2e} {:16.4f} {:14.4f}\".format(M, sample_bias_theta, sample_var_theta, true_var_theta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Uninformative prior and Laplace's principle of indifference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\theta$ is a probability, so it is bounded and must belong to the interval $[0, 1]$. We could assume that all values of $\\theta$ in $[0, 1]$ are equally likely. Thus our prior could be that $\\theta$ is uniformly distributed on $[0, 1]$, i.e. $\\theta \\sim \\mathcal{U}(a = 0, b = 1)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This assumption would constitute an application of **Laplace's principle of indifference**, also known as **principle of insufficient reason**: when faced with multiple possibilities, whose probabilities are unknown, assume that the probabilities of all possibilities are equal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, this principle should be used with great care, as we are assuming something strictly greater than we know. Saying \"the probabilities of the outcomes are equally likely\" contains strictly more information that \"I don't know what the probabilities of the outcomes are\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If someone tosses a coin and then covers it with her hand, asking you, \"heads or tails?\", it is probably relatively sensible to assume that the two possibilities are equally likely, effectively assuming that the coin is unbiased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If someone asks you, \"Is So-and-So a murderer?\", you should think twice before applying Laplace's principle of indifference and saying \"Well, it's a 50% chance that So-and-So is a murderer, it may be safer to lock So-and-So up.\" (Poor So-and-So!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the context of Bayesian estimation, we may be OK to apply Laplace's principle of indifference. This constitutes what is known as a **uninformative prior**. Our goal is, however, not to stick with a prior, but use the likelihood to proceed to a posterior based on new information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pdf of the uniform distribution, $\\mathcal{U}(a, b)$, is given by\n",
    "$$p(\\theta) = \\frac{1}{b - a}$$\n",
    "if $\\theta \\in [a, b]$ and zero elsewhere. In our case, $a = 0, b = 1$, and so, assuming $\\theta \\in [0, 1]$, our uninformative uniform prior is given by\n",
    "$$p(\\theta) = 1.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us derive the posterior based on this prior assumption. Bayes's theorem tells us that\n",
    "$$\\text{posterior} \\propto \\text{likelihood} \\cdot \\text{prior}.$$ In terms of our pdfs, this is\n",
    "$$p(\\theta \\, | \\, x_{1:n}) \\propto p(x_{1:n} \\, | \\, \\theta) p(\\theta) = \\theta^{\\sum x_i} (1 - \\theta)^{n - \\sum x_i} \\cdot 1.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $\\propto$ stands for \"proportional to\", so we may be missing a normalising constant. However, by looking at the shape of the resulting pdf (note, the function's argument is now $\\theta$, not $x_i$, so it is not the pdf of a Bernoulli distribution!), we recognise it as the pdf of the distribution\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\sum x_i, n - \\sum x_i\\right),$$\n",
    "and we immediately know that the missing normalising constant factor is\n",
    "$$\\frac{1}{B\\left(\\sum x_i, n - \\sum x_i\\right)} = \\frac{\\Gamma\\left(\\sum x_i\\right) \\Gamma\\left(n - \\sum x_i\\right)}{\\Gamma(n)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now assume that we have tossed the coin ten times and, out of those ten coin tosses, we get heads on seven. Then our posterior distribution becomes\n",
    "$$\\theta \\, | \\, x_{1:n} \\sim \\text{Beta}(\\theta \\, | \\, 7, 3).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, from the properties of this distribution,\n",
    "$$\\mathbb{E}[\\theta \\, | \\, x_{1:n}] = \\frac{\\sum x_i}{\\sum x_i + (n - \\sum x_i)} = \\frac{\\sum x_i}{n} = \\frac{7}{7 + 3} = \\frac{7}{10} = 0.7,$$\n",
    "$$\\text{Var}[\\theta \\, | \\, x_{1:n}] = \\frac{\\left( \\sum x_i \\right) \\left( n - \\sum x_i \\right)}{\\left( \\sum x_i + n - \\sum x_i \\right)^2 \\left( \\sum x_i + n - \\sum x_i + 1 \\right)} = \\frac{n \\sum x_i - \\left( \\sum x_i \\right)^2}{n^2 (n + 1)} = \\frac{7 \\cdot 3}{(7 + 3)^2 (7 + 3 + 1)} = \\frac{21}{1100} \\approx 0.019,$$\n",
    "and the posterior pdf looks as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "posterior mean: 0.7\n",
      "posterior s.d.: 0.13816985594155148\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVxV1fo/8M9iFgSRyYEDAgriAIKg4oxpV3MqU9OcIjOnvKbdurf63m+/7q3vvd0sK9NQ0zTNcig1NdOSNDVxQEEQFVQEPYggk8zjWb8/GC6w9+Ec4Jyz9zk879eLl7DWPmc/GxcPm73XfhbjnIMQQojxM5M6AEIIIbpBCZ0QQkwEJXRCCDERlNAJIcREUEInhBATYSHVjl1cXLiXl5dUuycSSMpJAgD0du4tcSTykfKoGADg42oncSTyQmNFvcuXL2dzzl3F+iRL6F5eXoiJiZFq90QC4dvDAQCnIk5JGoeczNoUDQDYs2SoxJHIC40V9Rhjaer66JILIYSYCMnO0En78/dRf5c6BGIkaKy0DiV0YjDjfMZJHQIxEjRWWkdWCb2yshJKpRJlZWVSh0L0oKK6AgBgZW4lWQw2NjZQKBSwtLSULAaiWdzDOABAUNcgiSMxLrJK6EqlEvb29vDy8gJjTOpwiI4lZdfOXHCRZuYC5xw5OTlQKpXw9vaWJAainVXHVgGgm6ItJauEXlZWRsmc6A1jDM7Oznj06JHUoegd5xyPyx+jsLwQVaoqVKmqwMHR0aoj7K3sYWdlBzNGcyJMjawSOgBK5kSvTHF8VamqcO/xPdzOvY07uXeQUZSBnJIcVKmq1L6GMYauHbvCs5MnenTqAV9nX3g4eJjk96c9kV1CJ4RoVl5VjmtZ13Al4woSshJQXlXeotdzzpFRmIGMwgxcUF4AADh1cEJQ1yCEdA9Bz849KbkbIUrorfTOO+9g1KhRGDeO7sYTw1EWKHEq9RQuKC/U32TWldzSXPx29zf8dvc3KBwUeML7CQx2HwxLc7qBbCwoobdCdXU1/vnPf7b4Nebm5nqKyDi4O7hLHYJR4pwj8VEijt0+hls5twyyT2WBEjuu7sCBmwcwyXcSRvUYBXMzw43ff439l8H2ZUpkndAPHwaOHNHNe23apN12qampmDBhAoYMGYLY2Fj4+flhx44d6Nu3LxYuXIhffvkFK1aswLFjxzB58mTMmDEDUVFReP3111FVVYVBgwYhMjIS1tbW8PLyavSa2bNn6+ZgjFRHq45Sh2BUOOdIyknCjzd/REpeSotfb2NhA6cOTrAyt4KFmQU4OIoqilBYXoiSyhKt3qOwvBC7r+3Gb3d/w/S+0zGgywCDXIoZ5jFM7/swRbJO6FJJSkrC1q1bMXz4cCxcuBBffPEFgJo5zGfPngUAHDt2DEDNzJyIiAhERUXBz88PCxYsQGRkJFatWiV4TXtXVFEEgBK7NjKLMrH72m5cf3Rdq+1tLW3h5+yHnk490bNzT3Tt2BW2lrZqk29xRTHuPb6HtMdpuJl9E0nZSVBxldr3zyrOQuSlSAR0CcD8wPnoZNOpVcelrXP3zwGgxN5SlNBFeHh4YPjw4QCAefPmYd26dQCAWbNmCbZNSkqCt7c3/Pz8AAAvvPACNmzYUJ/QxV7TXqUXpAOQbh66MSivKsfRW0fxa8qvqFZVN7utjYUNQruHYmC3gfB38W/RJRE7Kzv0ce2DPq59MKHXBJRUliA+Mx6n007jTu4dta9LyEzAu6fexZyAORjkPkjr/bXU21FvA6B56C1FCV1E07Oauq/t7IQlTjUtsi32GkLEJOckY3vcduSU5DS7XXf77hjjPQZD3IfA2sJaJ/u2tbRFmCIMYYowpOan4kTKCVxKvyS6bUllCbZc2YKErATMC5wn6ZO/pDFK6CLu3buH6OhoDB06FN999x1GjBiB2NhY0W39/f2RmpqK27dvo1evXti5cydGjx5t4IiJMePg2Ju4F1EpUc1up3BQYGrvqQjsEqjX69hejl5YNHAR/tTzT9iXuA/JOcmi211QXsCDwgdYFroMzrbOeouHaE/WCX3KlJoPQ+vTpw++/vprLFmyBL6+vli2bBk+//xz0W1tbGywbds2zJw5s/6m6NKlSw0cMTFWFdUVeFD4AFEpp9Ru42zrjGf7PIuQbiEGnRvu2ckTrw19DXEP4/BtwrcoKC8QbHP/8X3868y/sCR0Cfyc/QwWGxEn64QuFTMzM2zcuLFRW2pqaqOvt2/fXv/52LFjRc/gm76GkIbO3T+HtMdpai/bWZhZYEKvCZjQa4Jkc8EZYwjuFgxfZ198m/AtLj+4LNimqKIIn53/DIsGLkJwt2AJoiR1KKETg/Ho5CF1CLJQWV2JbxO+xbn758C5+CwOX2dfLBiwAG52bgaOTlxHq45YHLIY57ucx86rOwVlBapUVdh0eRMWDFigk5kpn074tM3v0R5RQm/Cy8sL165dkzoMk2RraSt1CJIrKC/AF5e+wN28u6L9FmYWmNZnGp7wfkKWxbPCFGHo1rEbImMikVea16iPc46v475GaWUpxvqMbdN+qGxu68hvxBCTVVBeIHodtr2ou96sLpl3t++O/xn1PxjnM06WybxOD8ceeHvk2+jp1FO0f2/iXpy8e7JN+ziRcgInUk606T3aIzpDJwaTUZgBAHCwdpA4EsOLz4zHl5e/VFt/JUwRhjkBc3Q2DVHfHKwdsCpsFTbGbERiVqKgf/e13bCxsMFQj9Ytfv3+6fcB0MpFLSXf0wBCTMQf9/5A5KVI0WTOGEOXjl0QERRhNMm8jpW5FZYPWo7Q7qGi/V9f/RqxGeLTfYl+UEInRE845zh66yh2XN0h+li9vbU9PBwU6GTdyWhL1VqYWeClgS+JnolzzvHllS8NVlCMUELXi4MHD+L6de1qcDR06NAhfPDBB3qIiBga5xz7ru/Djzd/FO13d3DHWyPego1FBwNHpntmzAwLBizAwG4DBX3VqmpExkTiUbHprxIlB5TQ9aA1Cb2qqgpTp07Fm2++2aLXEPnhnOPbhG/VPvnZx7UP/jr8ryb1dKUZM8NLA19CP7d+gr7iimJsuLRB6wqPpPVkeVN0yeElet/Hpini9XTVlc+Njo4WLZH75ptv4tChQ7CwsMCf/vQnPPvsszh06BB+//13vP/++/jhhx8AAK+88goePXoEW1tbfPnll/D390dERAScnJwQGxuLgQMHIiAgADExMVi/fj3S0tKwcOFCPHr0CK6urti2bRs8PT0Fr/n444/1/r3SlR6dekgdgt6puApfx32N88rzov2D3QfjhaAXYGEmyx+9NrEws8DS0KX4JPoTQbnfjMIMbL68GSuHrNRqBs+myVrWuyaN0Bm6iKSkJCxevBjx8fFwcHDA2rVrERERgT179iAhIQFVVVWIjIxEbm4uDhw4gMTERMTHx+Pvf/87hg0bhqlTp2LNmjWIi4tDz549sXjxYnz++ee4fPkyPvroIyxfvrx+X8nJyThx4oQgMa9YsQILFixAfHw85s6di5UrV2p8jdzZWNrAxtJG6jD0RsVV2B63XW0yH+czDguDF5pkMq9Td6NU7K+PG49u4ODNg1q9T2+X3lSVsxUooYtoWj43KipKUCL39OnTcHBwgI2NDRYtWoT9+/fD1lb44ExRURHOnTuHmTNnIigoCEuWLEFGRkZ9/8yZM0VXMoqOjsacOXMAAPPnz29UU13da+Quvywf+WX5UoehF5xz7Li6o359zqam9J6CGX1nGO3Nz5awt7bHisErYGMh/OV9/PZxxGfGa3yPw0mHcTjpsD7CM2kaEzpjzIMxdpIxdoMxlsgYe1Vkm3DG2GPGWFztxzv6CdcwtP2hs7CwwMWLFzF9+nQcPHgQEyZMEGyjUqng6OiIuLi4+o8bN27U92tbXrdhTMZakjezKBOZRZlSh6FznHN8E/8Nou9Hi/Y/2+dZTPab3C6SeZ3u9t3xcsjLose8LXabxhLBH0d/jI+jjesvUDnQ5gy9CsBfOOd9AIQBeIUx1ldkuzOc86Daj5YtuCkzdeVzAeC7777DuHHj6kvkAqgvkVtUVITHjx9j4sSJ+PTTTxEXFwcAsLe3R2FhIQDAwcEB3t7e2LdvH4CaH/6rV69qjGHYsGHYvXs3AGDXrl0YMWKEzo+TtB3nHHsS9+DsPfFVqWb3n43xvcYbOCp56O/WH8/4PyNoL6kswebLmwX1YEjbabyYxznPAJBR+3khY+wGAHcALZ+XpyV1NywNpWn53M8++wxhYWGCErm5ubl4+umnUVZWBs45PvnkEwDA7Nmz8fLLL2PdunX4/vvvsWvXLixbtgzvv/8+KisrMXv2bAwYMKDZGNatW4eFCxdizZo19TdFifwcTj6s9jH35/o9hzHeYwwckbyM7zket3Ju4VpW4/pIqfmpOHjzIGb0nSFRZKaJaVpxp9HGjHkBOA2gP+e8oEF7OIAfACgBPADwOudc8DwwY2wxgMUA4OnpGZKWltao/8aNG+jTp09Lj0GnUlNTMXnyZCrQpQdJ2UkApF+CTlfj7ETKCexL3CfaN73vdPyp5580vsesTTV/Ce5Z0rpH5I1BcUUx3j/9PnJLcxu1M8bw2tDXROuoh28PB0BL0IlhjF3mnIs+nqv1TVHGWEfUJO1VDZN5rSsAenDOBwD4HIDorWzO+WbOeSjnPNTV1VXbXRMiO9H3o9Um82f8n9EqmbcXdlZ2WBK6RLDmKecc22K3obSyVKLITI9WCZ0xZomaZL6Lc76/aT/nvIBzXlT7+VEAlowxF51GaiBUPld/vB294e3oLXUYbZaYlYgdV3eI9o3vNR5P+T5l4Ijkz8vRC0/3flrQnluaiz2JewTtO6ftxM5pOw0RmknRZpYLA7AVwA3O+Vo123St3Q6MscG179v8bWzS7lhZWMHKwrgXFE7LT8Omy5tEa7OM7DES0/ynSRCVcXiy55PwdfYVtEffjxYU8fLo5EELorSCNmfowwHMB/BEg2mJExljSxljdYtnzgBwjTF2FcA6ALN5Sy7Ok3YhtzRXcB3VmGSXZOPzi5+jvKpc0BfSPQRzAua0q6mJLWXGzNRWlfw24dtGpQH2XNuDPdeEZ+6kedrMcjkLoNlRyjlfD2C9roIipqmuQJNTByeJI2m5ksoSrLuwDoXlhYK+3i69sTB4oawXpZALF1sXzOo3S3DJqqC8APtv7Me8wHkAgMiYSADArP6zDB6jMaMRSIgGVaoqbIzZKPpQlLuDO5aFLjPpx/l1bZjHMAR2CRS0n0k7Q6V224gSuh5Q+VzTUVc5sW7KZUOdO3TGyiEr0cHS+EvgGhJjTO3qTDvjd6KyulKCqEwDJXQ9oPK5puOXO7/gj3t/CNptLGywcshKONo4ShCV8evcobPoDeTMokz8fPtnCSIyDZTQG0hNTYW/vz9eeOEFBAYGYsaMGSgpqblRExUVheDgYAQEBGDhwoUoL6+5Mfbmm2+ib9++CAwMxOuvv45z587h0KFDeOONNxAUFIQ7d+7gzp07mDBhAkJCQjBy5EjcvHkTABAREYHXXnsNY8aMwd/+9jds374dK1asAACkpaVh7NixCAwMxNixY3Hv3j3R1zS0fft2PPPMM5gyZQq8vb2xfv16rF27FsHBwQgLC0Nubs0NSXXxHD58GEOGDEFwcDDGjRuHzMyaSwzvvvsuFi5ciPDwcPj4+GDdunV6/p+Qh6sPr2L/DcEsXZgxMywNXYru9t0liMp0jPYaDZ/OPoL247ePo1JFZ+mtIdsLf/84nIjrD3S7Qnzf7g74f1OEBfgbSkpKwtatWzF8+HAsXLgQX3zxBVasWIGIiAhERUXBz88PCxYsQGRkJBYsWIADBw7g5s2bYIwhPz8fjo6OmDp1KiZPnowZM2oeax47diw2btwIX19fXLhwAcuXL8dvv/0G4L+lcM3NzbF9+/b6OOrK577wwgv46quvsHLlShw8eFDwmqauXbuG2NhYlJWVoVevXvjPf/6D2NhYrF69Gjt27MCqVauwePFi0XhGjBiB8+fPgzGGLVu24MMPP6wv0Xvz5k2cPHkShYWF6N27N5YtWwZLS8sWff/FfnjlKr0gHVtjt4r2zQmYgz6u0j7RbArMmBnmD5iP90+/j2pVdX17laoKT/d+GguDF0oYnXGSbUKXStPSuevWrcOTTz4pKJ+7YcMGrFixor587qRJkzB58mTB+zUsn1un7uweaL587v79NWeH8+fPx1//+leNrwGAMWPGwN7eHvb29ujUqROmTJkCAAgICEB8fHyz8SiVSsyaNQsZGRmoqKiAt/d/HwKaNGkSrK2tYW1tDTc3N2RmZkKhUGj4bjZmad6yXwBSKSwvxIZLG0SnJ47zGYeRPUZKEJVp6m7fHU/6PIljt481ar+TeweZRZlwsTXK5xMlI9uErulMWl+aziNmjEHdlPq68rlRUVHYvXs31q9fX3/mXadh+Vwxui6fa2393xtNZmZm9V+bmZmhqqqq2Xj+/Oc/47XXXsPUqVNx6tQpvPvuu6Lva25u3qrr99kl2QAg6x/SalU1Nl/eLFreNaBLAKb3nS5BVKZtou9EnFeeb1QrPyknCW9FvYX9s/bTDKIWoGvoTTQtnTtixAj4+/ubTPnc5uJ5/Pgx3N3dAQBff/21TvbXUE5JjsY62FLbd30fknOSBe3d7Lth0cBFNNdcD6wtrAVVF5NzkhHzIEbtuqxEHI3OJupK5wYGBiI3NxfLli2DjY0Ntm3bhpkzZyIgIABmZmZYunQpCgsLMXnyZAQGBmL06NGNyueuWbMGwcHBuHPnDnbt2oWtW7diwIAB6NevH378UXwl+IbWrVuHbdu2ITAwEDt37sRnn32ms2NUF8+7776LmTNnYuTIkXBxke9ZtL6cu39OtBSunZUdXhn0iugKPEQ3QruHilZdPHrrqOjDXERci8rn6lJoaCiPiYlp1CZ1+Vwqnatfci6fm5qfijV/rBEsumDGzLAqbJXeYm4P5XO1pSxQ4v9O/x9UXIXDyTXLz03xm4Jwr3A8H/C8xNHJh07K5xJiqgrLC7ExZqPoCjoz+82U/BdQe6FwUGBUj1GC9tNpp01y6UJ9oITeAJXObX9UXIUtV7YgrzRP0BemCMMYr/a94pChTfabLLi0peIqHLh5QKKIjAvdPiYG08upl9QhCPx480fczL4paO/h2APzAudR9UQDs7e2x/he41FUUdSoPTYjFrdybomW3yX/RWfoxGDMzcwFq9ZIKe5hnGD+MwB0tOqIpaFLjWbevKkZ5zMOLrYugumK+2/sVzuFmNSghE4MJqs4C1nFWVKHAaAmlm2xwoW3GWNYNHCRUZb4NRVW5lao5tVIfNR4WeKUvBQkZCVIFJVxoIRODCavNE/0WrWhcc6xKWYTyqrKBH1P936aHuuXgZgHMUgvSBe0/3jzRzpLbwYldD2g8rlCEyZMwKCeg7BkzhKdvWdERAS+//77Fr+uuLIYygKloH1A1wGY0GuCLkIjOiD2RLGyQInLGZcliMY4UELXAyqfK/TGG2/gP1/8R+owUFRRJFqjxdXOFRFBEXQTVEbsrOzg3Vm4qPihpEOia7oSSuiNUPlc/ZXPHTt2LOw6Nl+3Zt26dfXfy9mzZwv6OedYsWIF+vbti0mTJiErq2XX4yuqK0TXNLUws8CSkCWwtbRt0fsR/XvG/xlBW2ZRJi4oL0gQjfzJetpi+PZwQdtz/Z7D8kHLUVJZgom7Jgr6I4IiEBEUgeySbMzY27g+xKmIUxr3SeVz9Vc+V5MPPvgAd+/ehbW1NfLz8wX9Bw4cQFJSEhISEpCZmYm+ffti4ULtSqyquArZJdmi119n959NK8zLVG/n3vBz9hPU1zmSfASD3QfLataUHNAZehNNy+eePXsWSUlJgvK5p0+fhoODQ3353P3798PWVniG17BcbVBQEJYsWYKMjIz6/ubK586ZMwdATfncs2fPanwN8N/yua6uroLyuampqc3Go1QqMX78eAQEBGDNmjVITPzvLIO68rkuLi715XNbyrOTJzpadVTbHxgYiLlz5+Kbb76BhYXwXOP06dN4/vnnYW5uju7du+OJJ57Qet+5pbmiS5uFKcIwwlM3hc+I7pyKOIVTEafAGBM9S88uycaFdDpLb0rWZ+jNnVHbWto22+9i66LVGXlTVD63deVzDxw4gH/84x8AgC1btiA0VLTURLN++uknnD59GocOHcJ7772HxMREQWJvzTXuoooiFFcUC9q72XfDnIA5dN1c5no69UR/t/64ltX4Ke6fb/2MMEUYVcBsgL4TTVD53NaVz502bRri4uIQFxenNpnnluaK3pAEan7x3b9/H2PGjMGHH36I/Px8FBU1flpw1KhR2L17N6qrq5GRkYGTJ4WVEZtSd93c0twSi0MWiy5UTKT30bmP8NG5j+q/nuQ3SbBNVnEWYh7ECNrbM0roTVD5XP2Uzx05ciQWzV+EM7+fgUKhwPHjxxv1V1dXY968eQgICEBwcDBWr14NR0dHxMTEYNGiRQBqfmn4+voiICAAy5Ytw+jRo+tf/8477+DQoUON3rO56+ZzAubQmqAydiT5CI4kH6n/2qezj+jzAUdvHaV56Q1Q+dwGqHyufhm6fG5OSY6gJggAZKVlYdCAQQaJQRMqnyuubkJEw8umt3JuNTprr7M4ZDFCuocYKDLpUflc0u4UVxSLJnMLMwvYWWp334LIi6+zr9pFMOgsvQYl9AaofK5pqKyuRE6p+FJ3rnaudBPUiIldS1cWKAV1X9orSujEYBhjek+mnHO1182dOjjBytxKr/snutHBsgM6WHYQtPd27o2eTj0F7WJVM9sjSujEYPyc/UT/ZNal/LJ8VFRXCNptLW2bnQNP5OXnuT/j57k/C9oZY6L1dm7l3EJKXoohQpM1jQmdMebBGDvJGLvBGEtkjL0qsg1jjK1jjN1mjMUzxgbqJ1xC1CutLEVBeYGg3dzMHE4dnOhSi4kIcAsQnaFEZ+nanaFXAfgL57wPgDAArzDG+jbZ5ikAvrUfiwFE6jRKYhIeFD7Ag8IHennvalW12uvmLrYu9Ii4kXnv9/fw3u/vifYxxjC+13hB+9WHV5FRmCHyivZDY0LnnGdwzq/Ufl4I4AYA9yabPQ1gB69xHoAjY6ybzqM1ElQ+V1xheSEKy2seuLp58yaGDh0Ka2trfPSRcCpaS3DOkVOag2pVNYYHDkduzn8fJHK0cRSsUUnkL+puFKLuRqntH9R9EDp36Cxo/+XOL/oMS/ZadA2dMeYFIBhA0yIK7gDuN/haCWHSB2NsMWMshjEW8+jRo5ZFakTae/nc1NRUhIeHN7uNk5MT1q1bh9dff73N+yuqKEJpZamg3drCGg7WDm1+fyI/5mbmeNLnSUH7hfQLyC8TFnZrL7RO6IyxjgB+ALCKc970QqXYxUnBNAPO+WbOeSjnPNTV1bVlkRoAlc/VX/ncptzc3DBo0KBmKzZWV1cjIiIC/fv3R0BAQP2TuA1lZGVg0lOTMHHURLy16q36UWfGzOBi60LXzU3YCM8RgpLH1apqnEo9JU1AMqBVcS7GmCVqkvkuzvl+kU2UABrWH1UAaPPFUiqf277K5zYVFxeH9PT0+mcDmpbUVXEV/v7O3xEaFopX//oqfjv+G777+jsAgLOts2CRYWJarC2sMcZ7DH5K/qlR++m005joO7FdTlHVZpYLA7AVwA3O+Vo1mx0CsKB2tksYgMecc6O8O0Hlc1tfPnfatGkICgrCxIkTERMTg6CgIAQFBWHbtprFmM3NzFt0c9LHxwcpKSn485//jGPHjsHBofHlk/yyfET/EY1pz00DADwx/gl0cuwEO0s7WqzCyDnbOsPZ1lnjduFe4YJf3MUVxTivPK+v0GRNm1OY4QDmA0hgjNXVXH0bgCcAcM43AjgKYCKA2wBKALyoi+CofK54XHIsnwvUlNAFai5dRURE4NSpU436ezn1av4gm+jcuTOuXr2K48ePY8OGDdi7dy+++uorADVTFOtusDb83jAwOHZwbNF+iPz88NwPWm3nYO2Awe6Dce7+uUbtJ1JOYKTnyHZ3yU2bWS5nOeeMcx7IOQ+q/TjKOd9Ym8xRO7vlFc55T855AOfcaGtaUvnc1pXP1Yfs7GyoVCpMnz4d7733Hq5cuQKg5jppdkk2AGDIsCE4uK/mUtTJX08iPz+f6mO3M+N8xgnaMosykZCVIEE00qKR3wSVz9VP+VygpuaGskAJAHj48CEUCgXWrl2L999/HwqFAgUFje+1p6enIzw8HEFBQYiIiMC///1vcM7x8ecfY8fWHQCAV//2Ki6eu4hJoyfh0plL8PT0rH/9xIkT8eCBfua9E/1668RbeOvEW1pt6+7gLlpaNypF/bRHU0Xlcxug8rn6pYvyuQXlBcgrzRO021jYwM3OTas/saUeZw1R+VxxYuVzm3Mt6xo+v/C5oP1/R/8vFA4KHUYmPSqfS0xCRXWFaDKnKYqkn2s/dLMXPst48q7mVa1MCSX0Bqh8rnzVrT4kxtnWmR7tb+cYY3jCW7ho+IX0C6LryZoq2SV0KlRPxOSX5aOyulLQ3tGqY4umKNL4Ml1D3IcIxkJldSXO3jur5hWmR1YJ3cbGBjk5OfRDZ6IszS1had7yh5EaTlFs+n5i9TzU4ZwjJycHNjZU20XuFA6KFl/7trawxjCPYYL239N+h4qrdBWarMnqUTqFQgGlUglTrvNCgBsPb2i9rYqr8LjsseAHkjEGB2sHPDZ73KJ929jYQKEwrZtkpuibZ79p1evGeI9B1N2oRieFOSU5iM+MR1DXIF2FJ1uySuiWlpbw9vaWOgwiE5xzRMZE4upD4bz9mf1mYrjPcAmiInLmYuuCALcAxGfGN2o/efdku0josrrkQkzbqmOrsOrYKq23P3vvrGgy7+PaB2O9x+oyNCIzLR0rDY3xHiNou5l9U2+1+OWEEjoxmLiHcYh7KF4CoanMokzsTdwraLezskNEUARNUTRxLRkrTfVx6YMuHbsI2n9P/b2tYckeJXQiO1WqKmyN3Sq6Nuj8wPlwtKFaLUQ9xhjGeAnP0s8rz6O8qlyCiAyHEjqRnSPJR5CWnyZoH+E5AsHdgsDqx5wAAB0JSURBVCWIiBibMEWYoHxuWVUZLqZflCgiw6CETmQlOSdZdLFfNzs3PNfvOQkiIsaog2UHDFEMEbSfSj1l0tOiKaETg/Fz9oOfs5/a/pLKEnwV+5XgB86MmeGlgS/B2sJazSuJqdE0VrQxusdoQZuyQIm7+Xfb9L5yJqtpi8S0bZ6yWW0f5xzfJnwrWqtlau+p8HL00mNkRG6aGyva8ujkAZ/OPkjJS2nU/nvq7/Dp7NPm95cjOkMnsnAx/SIupV8StPs6+2J8r/ESRERMwWgv4Vl6zIMYFFUUSRCN/lFCJwaz+PBiLD68WNCeXZKNbxO+FbTbWNjgxaAXacGKdkjdWGmpkG4hsLNqvMJXlapKsMKRqaCfFGIwyTnJSM5JbtSm4ip8FfsVyqrKBNvPCZij1bqSxPSIjZXWsDS3FK3vcvbeWZO8OUoJnUjq51s/407uHUH7EMUQ0VkKhLTUSM+RgrbMokzcyr0lQTT6RQmdSCYlLwVHko8I2p1tnfF8/+cliIiYoi4du4iuknUm7YwE0egXJXQiibKqMmy9slW0iuKLQS+ig2UHiSIjpkjsLP1KxhWTW/yCpi0Sg2lY7e67hO9EVyB6qtdT8HX2NWRYRIZ0XRkxuFswOlp1bDS7pUpVhWhlNMb5jNPpvqRECZ0YzKcTPgUAXEq/hPPK84J+L0cvTPabbOiwiAzVjRVdsTCzwDCPYfjlzi+N2s+kncFY77EmU+yNLrkQg8opycGuhF2CdmsLa7w08CVaG5TozQjPEYK2h0UPcTv3tgTR6AcldGIwc/fPxZM7n0RpZamgb3b/2XCzc5MgKiJH8/bPw7z983T6nl06dhEtJ2BKa45SQicGE58Zj4dFDwXtId1DMFQxVIKIiFwpC5RQFih1/r4jewhvjl7OuCx6kmGMKKETg0jJS0FOSY6gvXOHzpgbMNdkrmESeQvuGgxbS9tGbZXVlbj0QFh2whhRQid6VzdFsSnGGBYGLxQ8mk2IvliaW4o+sGYql100JnTG2FeMsSzG2DU1/eGMsceMsbjaj3d0HyYxVpxz7IrfpXaKYltLpBLSUmI3R9Py0/RyicfQtDlD3w5ggoZtznDOg2o//tn2sIipuJB+oX6VmC52XdDFrmatR+/O3jRFkag1VDFUb/dVFA4K9HDsIWg3hbN0jfPQOeenGWNe+g+FmJqs4qxGVRQHuw8GUDtFMZimKBL1/j3u33p9/xGeIwTLHF5QXsD0PtNhaW6p133rk66uoQ9ljF1ljP3MGOunbiPG2GLGWAxjLObRo0c62jWRoypVFbZe2Sq6KO/cgLlwtXOVICpCagzqPkiQuEsqSxD7MFaiiHRDFwn9CoAenPMBAD4HcFDdhpzzzZzzUM55qKsr/UCbskNJh5Can9qo7ZeUXxCTEUNVFIlG0/dOx/S90/X2/h0sOyCkW4ig3djrpLc5oXPOCzjnRbWfHwVgyRhzaXNkxGhdf3Qdx28fF7SrVCrYWdKMFqJZTkmO6DRXXRK7OXoz+6be96tPbU7ojLGurHYSMWNscO17Gu93hLRJQXkBtsVuE7SbMTN0s+9Gqw8R2ejl1EvwdDLnHNHKaIkiajttpi1+ByAaQG/GmJIx9hJjbCljbGntJjMAXGOMXQWwDsBsbopLgRCNOOfYHrcdBeUFgr6pvafCxsJGgqgIEccYw1AP4Uyac/fPGe1qRhoTOuf8ec55N865JedcwTnfyjnfyDnfWNu/nnPej3M+gHMexjk37otQpNV+TfkViVmJgnZ/F39M6KVp5ishhjdUMVTwlHJOSY5Olr+TApXPJTqRkpeCAzcOCNo7WnXEwuCFYIxhrPdYCSIjxshQY6Vzh87o69pXcCJy7v450VWO5I4SOmmzksoSbLmyRbD6EAC8GPwiOtl0AgD87+j/NXRoxEgZcqwM8xgmSOiXMy5jdv/ZRrdyFt2hIm3COcfOqztFZwaM8xmH/m79JYiKEO0N6DJAtGBXzIMYiSJqPUropE1Op53GlYwrgnYvRy9M6zOtUdtTu57CU7ueMlRoxIgZcqxYmlvWP8XckDHOdqGETlrt/uP72Ju4V9BuY2GDl0NehoVZ4yt6pZWlJlN3muiXocfKMI9hgrY7uXeQWZRpsBh0gRI6aZWyqjJsvrwZVaoqQd/8AfPhYkvPlhHj4dnJE+4O7oJ2YztLp4ROWoxzjm/iv0FWcZagb2SPkQjtHipBVIS0HmNM9Cz9vPK86M1+uaKETlrs7L2zuJQuXOFF4aDArH6zJIiIkLYb7D5Y8CRzXmkekrKTJIqo5WjaImmR+4/vY/e13YJ2awtrLA5Z3GzpUap/TrQlxVhxsHZAf7f+iM+Mb9R+7v459HHtY/B4WoMSOtFaaWWp2uvmcwPmokvHLs2+/vVhr+srNGJipBorQz2GChJ67MNYlFaWGsWcdLrkQrTCOceOqztEr5sP9xxOJXGJSQjsEihY49aY5qRTQidaOZl6UnS+ucJBgef7P6/Ve4RvD0f49nAdR0ZMkVRjxcLMwqjnpFNCJxql5KVgX+I+QbuNhY3G6+aEGBuxtUzv5N4R/etUbiihk2YVlhdi8+XNolO35g+Yr/G6OSHGxrOTJ7rbdxe0n1eelyCalqGETtRScRW2xm5FXmmeoG+M9xiab05Mkro66eeV52VfJ50SOlHrcNJh3Hh0Q9Du3dkbM/rOkCAiQgxjiPsQ0Trpt3JvSRSRdmjaIhF19eFVHL11VNBuZ2WHxSGLBXVatPFcv+d0ERppB6QeK51sOonWSY++Hw0/Zz+JotKMEjoRyCrOwlexXwnaGWNYNHARnDo4tep9lw9a3tbQSDshh7EyVDFUbZ10awtriaJqHl1yIY2UV5Uj8lIkyqrKBH2T/Sajr2vfVr93SWUJSipL2hIeaSfkMFaCugYJHiYqrypH7MNYiSLSjBI6qcc5x874nXhQ+EDQF9AlAJN8J7Xp/SfumoiJuya26T1I+yCHsWJpbil64z/6vnznpFNCJ/VOpJwQLbrlYutSvy4oIe2J2Jz0pJwk0ZlfckAJnQAAbjy6gR9u/CBotzS3xLJBywRLdBHSHvh09oGbnVujNs65bOekU0InyC7JxpdXvhSdYzs/cD4UDgoJoiJEeowxhCnCBO1ynZNOCb2dq7sJWlxRLOgb5zOOim6Rdk8soT8seoi0x2kSRNM8mrbYjnHOsT1uO5QFSkGfv4s/pvedrtP9RQRF6PT9iOmS01hxtnWGn7MfknOSG7VH34+Gl6OXNEGpQQm9Hfvp1k+iFRSdOjhh0cBFgtVb2kpOP6RE3uQ2VsIUYYKEfunBJczsN7NVD9npC11yaafiHsbhcNJhQbuluSWWD1oOe2t7ne8zuyQb2SXZOn9fYnrkNlZCuocIqooWVxQjITNBoojEUUJvh5QFStEnQYGaMyOPTh562e+MvTMwYy/VgCGayW2s2FjYYGC3gYJ2udVJp4TezhSUF2D9xfUoryoX9D3l+xRVUCREDbGbowmZCSgsL5QgGnEaEzpj7CvGWBZj7JqafsYYW8cYu80Yi2eMCX+NEVmorK5E5KVI0YciArsE4uneT0sQFSHGwd/FH442jo3aVFyFSw+ED+NJRZur+dsBrAewQ03/UwB8az+GAIis/VdvDh8GjhzRzXtt2tR8/5IlutnP5MnAlCnq+/V9THWP9afkpeD0mcZ9dtXdUVX0EpZubdmToC09pmTnmn9b8z1tL/9PDbXnY9I0VqQ4JjNmhiGKITh++3ij/r+tj8bAoida9J6a4m8tjWfonPPTAHKb2eRpADt4jfMAHBlj3XQVINGNn279hAvKC4J2S94R/UpegQVsJIiKEOMiVgqgyPweis2E9Y+koIv5Nu4A7jf4WlnbltF0Q8bYYgCLAcDT01MHuybauJh+UXRGixnM0bd4GWxULgaJo2/xMoPshxg/uY6Vbvbd0MOxB9LyGz9UlGkVDZ8y3T630Rq6uCkq9ne66DOxnPPNnPNQznmoq6urDnZNNLmTewdfx30t2terZB46VfcyWCw9y2ahZ9ksg+2PGC85jxWxs/Qsq/PgEK67a2i6SOhKAA3nuSkAyOPvj3YusygTGy5tQJWqStDnWf4UulYOM2g8RWb3UWR2X/OGpN2T81gZ5D4I5mbmjdoqWAHyLK5LFNF/MW0KzDDGvAAc4Zz3F+mbBGAFgImouRm6jnM+WNN7hoaG8piYmJbGS7RUWF6ID85+IPpwRkj3ELw88GWDl8MN3x4OADgVccqg+5WzWZtq5jHvWSI862vP5D5WIi9FIu5hXKO20O6heDnkZb3vmzF2mXMuOr9Y4zV0xth3AMIBuDDGlAD+HwBLAOCcbwRwFDXJ/DaAEgAv6iZs0loV1RXYcGmDaDL37uyNF4NepNrmhLTBUI+hgoQe9zAOJZUlkpaa1pjQOefPa+jnAF7RWUSkTVRchS1XtuBu3l1Bn4utC5YPWi54hJkQ0jL93fqjo1VHFFUU1bdVqaoQ8yAGo3qMkiwuelLUhHDO8W3Ct7j68Kqgz87KDiuHrISDtYMEkRFiWizMLDDYXXhlWerl6Sihm5Cfbv2EM2lnBO0WZhZYPmg5unTsIkFUhJimYR7CSQUpeSnILMqUIJoa8qn7SNrkTNoZ0bnmAPBi8Ivo5WS46Ynq/GXoX6QOgRgJYxgrCgcFFA4KwXoC0cpoPOP/jCQx0Rm6CbiScQW7EnaJ9j3X7znZFNya0nsKpvTWw/POxOQYw1hhjGGoh3B2UvT9aKi4NHPSKaEbuRuPbmDrla2i6xuO7zUeY33GShCVuKTsJCRlJ0kdBjECxjJWBrsPFiwEk1+Wj5vZNyWJhxK6EbubdxeRMZGiDw6FKcIwzX+aBFGpt+TIEiw5oqOKU8SkGctYcbB2QH83weM5OHf/nATRUEI3WukF6fj84ueidc37u/XHggELaK45IQYgdnO0bk66oVFCN0KZRZn49PynKK4oFvT1cuqFJaFLBI8mE0L0I6BLAOys7Bq1VVZXIuaB4Z+Ep4RuZHJKcvDJ+U9QUF4g6FM4KPDK4FdgZW4lQWSEtE8WZhYY4i5cAuKPe38YPBZK6EYkrzQPa6PXiq445GbnhlfDXpX0sWNC2iuxyy6p+al4UGjYOoU0D91I5JflY230WtH6LE4dnLB66GrZPwX691F/lzoEYiSMbax4dPKARycP3H/cuELkH/f+wMx+Mw0WB52hG4GC8gKsjV6LrOIsQZ+DtQNWD10Npw5OEkTWMuN8xmGczzipwyBGwBjHynCP4YK2C+kXUK2qNlgMlNBlri6Ziz1ObGdlh9VDV8PNzk2CyFou7mGcoEIdIWKMcawMdh8MC7PGFz0KywuRkJVgsBgooctYflk+Pjr3ETIKBav5wdbSFqvDVqO7fXcJImudVcdWYdWxVVKHQYyAMY4VOys7BHUNErQbck46JXSZyi/Lx8fnPhY9M7e1tMXqoavh0clD5JWEEKmI3RxNyEzA47LHBtk/JXQZyinJwUfnPhK9Zm5jYYNXw16FZydaZJsQuenj2geONo6N2lRchWilYcrqUkKXmaziLKw5twaPih8J+urOzL0cvQwfGCFEIzNmJnqWfvbeWdF6Szrfv973QLSWXpCONX+sEZ1nbmtpi1VhqyiZEyJzwz2Fs10eFT/Crdxbet83zUOXibt5d/H5xc9FH+e3s7LD6jDjv2b+r7H/kjoEYiSMeay42Lqgj2sf3Hh0o1H7mbQz8HP20+u+KaHLwPVH17ExZqNooS0HawesClsFdwd3CSLTLbE/RQkRY+xjZYTnCEFCv5JxBbMrZgvqvugSXXKRWMyDGKy/uF40mXfu0BmvD3vdJJI5UDN9S6qyosS4GPtYCeoaJEjcVaoqXEy/qNf9UkKX0G93f8OWK1tEnyRztXPFG8PeMKl1QN+OehtvR70tdRjECBj7WLEws8BQhXA1ozP3zuj15igldAlwzvH99e+x59oe0f9chYMCfx3+VzjbOksQHSFEF0Z4jhC0pRek427+Xb3tkxK6gVWpqrDlyhb8eudX0X5fZ1/8ZdhfZF9oixDSvG723dDTqaeg/XTaab3tkxK6ARWWF2Jt9Fq1he8HdB2AV4dQCVxCTMWoHqMEbTEPYkRns+kCJXQDeVj0EB+c/QB3cu+I9o/qMQpLQ5fC0tzSwJERQvQlpFuI6GpG55Xn9bI/mrZoADce3cCmy5tQWlkq2v+M/zOY0GuCya8B+umET6UOgRgJUxkrluaWGKoYihMpJxq1/572O57wfkLnP/OU0PWIc47f7v6Gfdf3id78NDczx4IBCxCmCJMgOsMTq0RHiBhTGiujeowSJPTMokwk5ySjt0tvne6LLrnoSZWqCjuu7sDexL2iybzuUf72kswB4ETKCcHAJkSMKY2VLh27wN/FX9Cuj5ujWp2hM8YmAPgMgDmALZzzD5r0RwBYAyC9tmk953yLDuM0KnmledgYsxGp+ami/W52blgxeIVJzTHXxvun3wcAo1uJhhieqY2VUT1G4Wb2zUZtVzKuoKC8QKcz2jQmdMaYOYANAJ4EoARwiTF2iHN+vcmmezjnK3QWmZG6mX0TX17+EkUVRaL9/i7+WByyWK+P/xJC5CWoaxAcrB1QUF5Q36biKvxx7w885fuUzvajzSWXwQBuc85TOOcVAHYDeFpnEZgIzjl+vvUzPj3/qdpkPsZ7DFYOWUnJnJB2xtzMXPRBI12X1dXmkos7gIZLWSsBDBHZbjpjbBSAZACrOef3m27AGFsMYDEAeHqazgINheWF2Ba3DYlZiaL95mbmeL7/8xjZY6SBIyOEyMUIzxH4+fbP4Jyjk00njPAcgRGeI3Q600WbhC62t6a/Ug4D+I5zXs4YWwrgawBPCF7E+WYAmwEgNDRU/9XeDeBWzi1subIF+WX5ov2ONo5YGroU3p29DRwZIUROnG2dMdF3IhQOCgzoMgDmZuY634c2CV0JoGEhbgWABw034JznNPjySwD/aXto8latqsaR5CP1v3HF9HbpjZcHvgx7a3sDRydPmyZvkjoEYiRMdaxM7T1Vr++vTUK/BMCXMeaNmlksswHMabgBY6wb57xuafqpABoXAjYx2SXZ2HplK1LyUtRuM6HXBDzt/zTMGM0MraPrObfEdNFYaR2NCZ1zXsUYWwHgOGqmLX7FOU9kjP0TQAzn/BCAlYyxqQCqAOQCiNBjzJLhnOPsvbPYd32faP1yoGZ1oYXBC9Hfrb+Bo5O/w0mHAQBTek+ROBIidzRWWkereeic86MAjjZpe6fB528BeEu3ocnL47LH2HF1B65lXVO7TU+nnnh54Mvo3KGzASMzHh9HfwyAfkiJZjRWWoce/deAc45oZTT2Je5DSWWJ6DaMMUz2m4yJvhPpEgshRDKU0JuRW5qLb+K/UTsdEai5c/1S8EuidY8JIcSQKKGLUHEVolKicDj5sNpr5QAw3HM4nuv3HGwsbAwYHSGEiKOE3sTdvLvYlbAL9x8LnouqZ29tjwUDFiCwS6ABIyOEkOZRQq9VWF6I/Tf2a1xpfIhiCGb1m0WP77fCzmk7pQ6BGAkaK63T7hN6laoKv6f+jsPJh9UuQAHUPPE5N3AunZW3gUcnD80bEQIaK63VbhM65xzxmfH4/vr3yCrOUrsdYwzhXuF4xv8ZulbeRnuu7QEAzOo/S+JIiNzRWGmddpnQU/JSsP/GftzKudXsdh6dPDAvcB68HL0ME5iJi4yJBEA/pEQzGiut064S+oPCBzh48yCuPrza7Ha2lraY2nsqRnuNpnnlhBCj0S4SekZhBn669RNiHsQ0W3uYMYZhHsMwzX8aFdQihBgdk07o6QXp+Pn2zxoTOVBTDGhm35l0M4YQYrRMMqGn5KXg2O1jGi+tAEDXjl3xbJ9nEdglUKeF5gkhxNBMJqGruApXH17Frym/4k7uHY3bd+7QGVP8pmCox1C6Tm4g3z/3vdQhECNBY6V1jD6hl1SW4I97f+BU6ilkl2Rr3N7B2gETek3AqB6jYGluaYAISR0XWxepQyBGgsZK6xhlQuecI+1xGs6kncGF9AuorK7U+Bp7a/v6RG5lbmWAKElT2+O2AwAigiIkjYPIH42V1jGqhF5aWYqL6Rdx5t6ZZmutNOTUwQnje43HcI/hdEYuMfohJdqisdI6RpPQjyQfwfHbx1FRXaHV9goHBZ7s+SQGdR+kl8VYCSFEbowmoXew6KBVMu/n1g9P+jwJfxd/mrVCCGlXjCahhynCsP/GflSpqgR9HSw7YJjHMIR7hcPNzk2C6AghRHpGk9DtrOwwsNtAXEy/WN/m09kHI3uMREi3EFhbWEsYHSGESM9oEjoAjOwxEglZCQhThGGE5wgoHBRSh0Ra4Ojco5o3IgQ0VlrLqBK6r5MvPnzyQ5p2aKRsLW2lDoEYCRorrWNUj0gyxiiZG7EvLn2BLy59IXUYxAjQWGkdo0roxLjtTdyLvYl7pQ6DGAEaK61DCZ0QQkwEJXRCCDERlNAJIcREUEInhBATwTSt5KO3HTP2CECaJDtvHRcAmuvzyp+pHAdgOsdiKscBmM6xyPk4enDOXcU6JEvoxoYxFsM5D5U6jrYyleMATOdYTOU4ANM5FmM9DrrkQgghJoISOiGEmAhK6NrbLHUAOmIqxwGYzrGYynEApnMsRnkcdA2dEEJMBJ2hE0KIiaCETgghJoISehOMsQmMsSTG2G3G2JvNbDeDMcYZY7Kc2qTpOBhjEYyxR4yxuNqPRVLEqYk2/x+MsecYY9cZY4mMsW8NHaO2tPg/+aTB/0cyYyxfiji1ocWxeDLGTjLGYhlj8YyxiVLEqYkWx9GDMRZVewynGGPyXoSBc04ftR8AzAHcAeADwArAVQB9RbazB3AawHkAoVLH3ZrjABABYL3UsergOHwBxALoXPu1m9Rxt2VsNdj+zwC+kjruNvy/bAawrPbzvgBSpY67lcexD8ALtZ8/AWCn1HE390Fn6I0NBnCbc57COa8AsBvA0yLbvQfgQwBlhgyuBbQ9DrnT5jheBrCBc54HAJzzLAPHqK2W/p88D+A7g0TWctocCwfgUPt5JwAPDBiftrQ5jr4Aomo/PynSLyuU0BtzB3C/wdfK2rZ6jLFgAB6c8yOGDKyFNB5Hrem1f0p+zxjzMExoLaLNcfgB8GOM/cEYO88Ym2Cw6FpG2/8TMMZ6APAG8JsB4moNbY7lXQDzGGNKAEdR8xeH3GhzHFcBTK/9fBoAe8aYswFiaxVK6I0xkbb6eZ2MMTMAnwD4i8Eiap1mj6PWYQBenPNAACcAfK33qFpOm+OwQM1ll3DUnNVuYYw56jmu1tDmWOrMBvA957xaj/G0hTbH8jyA7ZxzBYCJAHbW/vzIiTbH8TqA0YyxWACjAaQDqNJ3YK0lt2+w1JQAGp6pKtD4T0V7AP0BnGKMpQIIA3BIhjdGNR0HOOc5nPPy2i+/BBBioNhaQuNx1G7zI+e8knN+F0ASahK83GhzLHVmQ76XWwDtjuUlAHsBgHMeDcAGNQWv5ESbn5MHnPNnOefBAP6ntu2x4UJsGUrojV0C4MsY82aMWaHmB+tQXSfn/DHn3IVz7sU590LNTdGpnPMYacJVq9njAADGWLcGX04FcMOA8WlL43EAOAhgDAAwxlxQcwkmxaBRakebYwFjrDeAzgCiDRxfS2hzLPcAjAUAxlgf1CT0RwaNUjNtfk5cGvxl8RaArwwcY4tQQm+Ac14FYAWA46hJcHs554mMsX8yxqZKG532tDyOlbXT/K4CWImaWS+youVxHAeQwxi7jpqbVm9wznOkiVi9Foyt5wHs5rXTKuRIy2P5C4CXa8fXdwAi5HZMWh5HOIAkxlgygC4A/k+SYLVEj/4TQoiJoDN0QggxEZTQCSHERFBCJ4QQE0EJnRBCTAQldEIIMRGU0AkhxERQQieEEBPx/wGVRCwdkMuwygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alpha, beta = 7, 3\n",
    "\n",
    "posterior_mean, posterior_var, posterior_skew, posterior_kurt = scipy.stats.beta.stats(alpha, beta, moments='mvsk')\n",
    "\n",
    "xs = np.linspace(scipy.stats.beta.ppf(0.01, alpha, beta), scipy.stats.beta.ppf(0.99, alpha, beta), 100)\n",
    "plt.plot(xs, [1 for x in xs], 'b--', lw=5, alpha=.6, label='prior')\n",
    "plt.plot(xs, scipy.stats.beta.pdf(xs, alpha, beta), 'g-', lw=5, alpha=.6, label='posterior')\n",
    "plt.axvline(posterior_mean, label='posterior mean')\n",
    "posterior_sd = np.sqrt(posterior_var)\n",
    "plt.axvline(posterior_mean - posterior_sd, linestyle='--', color='g', label='posterior mean - 1 s.d.')\n",
    "plt.axvline(posterior_mean + posterior_sd, linestyle='--', color='g', label='posterior mean + 1 s.d.')\n",
    "plt.legend(loc='upper left');\n",
    "\n",
    "print('posterior mean:', posterior_mean)\n",
    "print('posterior s.d.:', posterior_sd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the mean of the posterior, 0.7, matches the frequentist maximum likelihood estimate of $\\theta$, $\\hat{\\theta}_{\\text{ML}}$, and our intuition. Again, it is not unreasonable to assume that the probability of getting heads is 0.7 if we observe heads on seven out of ten coin tosses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A more informative prior: the Beta distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us question our prior. Is it somewhat *too* uninformative? After all, most coins in the world are (probably!) close to being unbiased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could use a $\\text{Beta}(\\alpha, \\beta)$ prior instead of the Uniform prior. Picking $\\alpha = \\beta = 2$, for example, will give a distribution on $[0, 1]$ centred on $\\frac{1}{2}$, incorporating a prior assumption that the coin is unbiased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pdf of this prior is given by\n",
    "$$p(\\theta) = \\frac{1}{B(\\alpha, \\beta)} \\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1}$$\n",
    "on $\\theta \\in [0, 1]$, and so the posterior becomes\n",
    "$$p(\\theta \\, | \\, x_{1:n}) \\propto p(x_{1:n} \\, | \\, \\theta) p(\\theta) = \\theta^{\\sum x_i} (1 - \\theta)^{n - \\sum x_i} \\cdot \\frac{1}{B(\\alpha, \\beta)} \\theta^{\\alpha - 1} (1 - \\theta)^{\\beta - 1} \\propto \\theta^{(\\alpha + \\sum x_i) - 1} (1 - \\theta)^{(\\beta + n - \\sum x_i) - 1},$$\n",
    "which we recognise as a pdf of the distribution\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha + \\sum x_i, \\beta + n - \\sum x_i\\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why did we pick this prior distribution? Its pdf lives entirely on the compact interval $[0, 1]$, unlike, for example, the normal distribution, which has tails extending to $-\\infty$ and $+\\infty$. With the parameters $\\alpha = \\beta = 2$, it is centered on $\\theta = \\frac{1}{2}$, incorporating the prior assumption that the coin is unbiased."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we initially assume a $\\text{Beta}(\\theta \\,|\\, \\alpha = 2, \\beta = 2)$ prior, then the posterior expectation is\n",
    "$$\\mathbb{E}\\left[\\theta \\,|\\, x_{1:n}\\right] = \\frac{\\alpha + \\sum x_i}{\\alpha + \\sum x_i + \\beta + n - \\sum x_i} = \\frac{\\alpha + \\sum x_i}{\\alpha + \\beta + n} = \\frac{2 + 7}{2 + 2 + 10} = \\frac{9}{14} \\approx 0.64.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that both the prior and posterior belong to the same probability distribution family &mdash; Beta. In Bayesian estimation theory we refer to such prior and posterior as **conjugate distributions** (with respect to this particular likelihood function)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unsurprisingly, since now our prior assumption is that the coin is unbiased, $\\frac{1}{2} < \\mathbb{E}\\left[\\theta \\,|\\, x_{1:n}\\right] < \\frac{7}{10}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps surprisingly, we are also somewhat more certain about the posterior (its variance is smaller) than when we assumed the uniform prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "posterior mean: 0.6428571428571429\n",
      "posterior s.d.: 0.12371791482634838\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeVzVVf7/n+fCBURRVNwBAfcFBEVTM8Wt/OZSmZYtU4yVZdlMMzXfmek30/6drakp252pbLHRyhbbzDSVLDcUcgVXQAQVkFV27vn98VHk8rnABe5+z/PxuA8u78/5fD5vj5f3PZ9z3uf1FlJKFAqFQuEZGJztgEKhUChshwrqCoVC4UGooK5QKBQehArqCoVC4UGooK5QKBQehK+zbhwSEiIjIiKcdXuFh5NekA7AkO5DnOyJ9ZzIuwBAVI+OTvbEetyxn92dPXv25EspezR13GlBPSIiguTkZGfdXuHhJKxMAGBL4han+tEabn5jOwBr7p3gZE+sxx372d0RQmQ2d1xNvygUCoUH4bSRukJhT/40+U/OdsErUP3seqigrvBIZkTNcLYLXoHqZ9fDpYJ6TU0N2dnZVFZWOtsVhZtTXVcNgJ+Pn5k9ICCA0NBQjEajM9zyOFLPpAIQ2zvWyZ4oLuFSQT07O5ugoCAiIiIQQjjbHYUbk55/MSsj5HJWhpSSgoICsrOziYyMdJZrHsVD6x8C1EKpK+FSQb2yslIFdIXdEELQvXt38vLynO1Kq6mpq+FM2RkqaysJ8A3A39efTn6dCDQGOts1hYvhUkEdUAFdYVfc5fMlpeTo+aPsOr2Lk4UnySnNwSRNunZhXcKI7hlNTK8YIoLVgEjhgkFdofBmJJIfs37k+5Pfk12S3WL7U8WnOFV8iq+Pfk2/zv2YP2w+I3qMUMHdi1F56m3kscceY+PGjc52Q+FBlNeUk1GUwbs/v2tVQG/M6ZLTvLTzJV7Y8QKnik/ZwUOFO6BG6m2grq6Op556qtXn+Pj42MkjRWP6de7nbBesprK2krWH1pJdUmuT66Xlp/GXH/7CzSNvZkr/KXYdtf9l+l/sdm1F23DpoP7FF/Dll7a51htvWNcuIyODWbNmccUVV5CSksLgwYN59913GT58OIsXL2bDhg0sW7aM9evXM2fOHBYsWMCmTZt45JFHqK2tZezYsbz22mv4+/sTERFhds6iRYts849RtEgnv07OdsEqckpzeHnXyxSUFwATm20bHBBMtw7dqK6rpqK24uI5ljFJE//d/18yizK5NfpWjD72SeGcGNa8zwrH02JQF0IEAEmA/8X2H0spH2/Uxh94FxgDFAA3SykzbO6tg0hPT+fNN9/kyiuvZPHixbz66quAluO8bds2ANavXw9oGTuJiYls2rSJwYMHc8cdd/Daa6/x0EMP6c5ROI6y6jLAtYN7Wn4arye/TkVNRZNtBnQbwNSIqQwJGUJn/85mx85XnOfAuQP8dOonThaetHj+T6d+Iqc0h2XjlhHkH2RT/y9dH1RwdyWsmVOvAqZJKUcBscAsIcT4Rm3uAgqllAOBfwF/t62bjiUsLIwrr7wSgNtvv70+KN988826tunp6URGRjJ48GAA7rzzTpKSkuqPWzpHYX9Ol5zmdMlpZ7vRJDuzd7J85/ImA/qg7oN49KpH+d8r/5ex/cbqAjpAtw7dmNx/Mr+/8vfcG38vPTv2tHitjKIMXtz5IuU15Tb9NwA8uulRHt30qM2vq2g7LQZ1qVF28VfjxVfjatXXAe9cfP8xMF248fJ7Y9cv/d6xo14StaXC3ZbOUXg3P2T+wFspb1FnqtMdMwgDi0Yu4uEJD9M/uL9V1xNCMLrPaB5PeJyEiASLbU4Vn+KlnS9RVVvVHtcVboBV2S9CCB8hRCpwDvhOSrmzUZN+wCkAKWUtUAx0t6WjjiQrK4vt2zUZ1P/+979MmjSpybZDhw4lIyODY8eOAfDee+8xZcoUh/ipcD+Sc5JZtX+VxWN+Pn70D+7P1MipbVrc9DX4ckv0Lfxi1C/wMegX5U8UnuDV3a9SU1fT6msr3AerFkqllHVArBAiGPhUCDFSSnmgQRNLn0DdEFYIsQRYAhAeHt7ifefO1V6OZtiwYbzzzjvce++9DBo0iKVLl/LSSy9ZbBsQEMDbb7/NwoUL6xdK77vvPgd7rHAHDp47yFspb1l8uhvYbSDnu4RhEO3PkJoUPoleHXuxfOfyeg2cS6Tlp/H+vvdJjE1UueweSquyX6SURUKILcAsoGFQzwbCgGwhhC/QBThv4fwVwAqA+Pj45uctnIjBYOD11183s2VkZJj9vnLlyvr306dPJyUlRXedxucovJeThSd5Lfk1i1Muo/uMZnHcYm4/YLuiMYO6D+L+sffz8q6XqTWZp0ruyN7BwG4Duar/VTa7n8J1aHH6RQjR4+IIHSFEB2AGkNao2TrgzovvFwDfy5YmmxUKOxLWJYywLmHOdgOAkqoSXk9+3eK0R0yvGO4efbddUg6H9RjGkjFLMAj9n/nqA6vJKs5q9z1emPUCL8x6od3XUdgOa+bU+wCbhRD7gN1oc+pfCiGeEkLMu9jmTaC7EOIY8FvgD/Zx1/5ERERw4MCBlhsqXJpAY6BLiF3VmepYsWcFRZVFumODug9iyZglFue/bcWo3qO4PeZ2nb3WVMsbyW+0OyMmtneskt11MVqcfpFS7gPiLNgfa/C+ElhoW9cUirZTUlUCYDEV0JF8dOgjjhYc1dnDuoTxwNgH7LYpqCFXhl/JicITbMsy3y+RX57P+/veZ8mYJW2+9sYTmlSGKpbhOijtF4VHkluaS25prlN92Jm9k80nN+vsQf5BPDD2AToYOzjMl0UjF1mcjtqTs4eUXP16kLU8k/QMzyQ90x7XFDZGBXWFwg4UlBfwwf4PdHaDMLBkzBK6dujqUH+MPkbuHXOvxS+S/x74r102JimcgwrqCoWNMUkTb6e+TWWtvizjguELGNx9sBO8gh4de3Bb9G06e3FlMWsPrXWCRwp7oIK6Hfjss884dOhQq89bt24df/vb3+zgkcKRbDyx0eI8+th+Y5kWOc0JHl0mvm88Mb1idPZtWdtIy2+c1KZwR1RQtwNtCeq1tbXMmzePP/zB+sSh2lrbSLUqbEd2STafpX2ms3ft0JVbo291+oYfIQS3Rt9KgG+A7tiqfat0Oe0K98MlpXfv/eJeu9/jjbmWtXibkt7dvn27RXndP/zhD6xbtw5fX1+uvvpq5s+fz7p169i6dSvPPPMMa9dqj7UPPPAAeXl5BAYG8u9//5uhQ4eSmJhIt27dSElJYfTo0URHR5OcnMzLL79MZmYmixcvJi8vjx49evD2228THh6uO+e5556ze1+5I/27WKebYktM0sTK1JW6DUZCCH4Z+0uXSLEE7QvmxuE3smqfuVzBuQvn2JqxlelR062+1htzrNS0VjgMNVK3QHp6OkuWLGHfvn107tyZ559/nsTERNasWcP+/fupra3ltdde4/z583z66accPHiQffv28ac//YmJEycyb948nn32WVJTUxkwYABLlizhpZdeYs+ePfzzn//k/vvvr7/XkSNH2Lhxoy44L1u2jDvuuIN9+/Zx22238atf/arFcxSXCTAGEGDUj0btycYTGy1WHJoRNYMhIUMc6ktLXBV+FYO6D9LZvzzyJReqL1h9nSEhQ1zu3+btqKBugcbSu5s2bbIor9u5c2cCAgK4++67+eSTTwgM1I/EysrK+Omnn1i4cCGxsbHce++95OZeTrVbuHChxYpI27dv59ZbbwXgF7/4hZkme1PnKC5TVFlkccOPvSgoL+CL9C909r5BfbluyHUO88NahBDcPOJm3XRQeU05Xx39yurrfJH+hcV/t8J5qKBuAWvnPX19fdm1axc33ngjn332GbNmzdK1MZlMBAcHk5qaWv86fPhw/XFrpXkb+qTkfFvmbNlZzpaddci9pJR8sP8DnXiWEII7Y+90yAajthDWJcxicYvNJzdz7sI5q67x3PbneG67emJ0JVRQt0Bj6d0ZM2ZYlNctKyujuLiYa6+9lhdeeIHU1FQAgoKCKC0tBaBz585ERkby0UcfAVoA+Pnnn1v0YeLEiaxevRqAVatWNSv/q3Aue3L3cOCcXloiISKBiOAIxzvUCuYNmYefj5+ZzSRNKsXRjXHJhdKmFjEdRWPp3RdffJHx48fr5HXPnz/PddddR2VlJVJK/vWvfwGwaNEi7rnnHpYvX87HH3/MqlWrWLp0Kc888ww1NTUsWrSIUaNGNevD8uXLWbx4Mc8++2z9QqnC9aiqreLDgx/q7MEBwVw/9HoneNQ6ggOCuWbgNboplNQzqWQWZVpdqEPhOghniSnGx8fL5GRzqdHDhw8zbNgwp/hziYyMDObMmaNEvdyc9Px0AIuLeLb8nH2W9hnfHP1GZ78v/j7i+ugkk5rl5je0p8M1906wiW/WUlVbxWObH9OtQcT0iuGBcQ80e27CygQAtiRusZN3isYIIfZIKeObOq6mXxSKNpJfns93x7/T2WN6xbiVcqG/rz9zBs/R2fed3UdmUaYTPFK0BxXUG6Gkdz2DyOBIIoMj7XqPtYfW6jbr+Bp8uXmkPqvE1ZkQNoHugfoKlF8e+bLZ89674T3eu+E9e7mlaAMqqCs8Ej9fP/x8/Vpu2EbS89PZm7tXZ585YCYhgSF2u6+98DX4cu2ga3X2lkbrrlSMRKGhgrrCIzlfcZ7zFbqKijbBJE0WF0e7BHRh1kB9Wqu7MD50fKtH62sOrGHNgTX2dEvRSlRQV3gkeRfyyLuQZ5dr7zq9i+ySbJ39hqE3WNRUcReaG61b+vcCvJb8Gq8lv2Zv1xStQAV1haIV1NTV8Hna5zp7RHAE40PHO8Ej29LUaN3SgrDCNVFB3Q4o6V3PZWvmVovTOguGL3C7xVFL+Bp8LU4h7Tq9i8KKQid4pGgtKqjbASW965lU1FTw9dGvdfaYXjEWxbHclfGh4+nk18nMZpImvj/5vZM8UrQGFdQbkJGRwdChQ7nzzjuJiYlhwYIFlJdrZb42bdpEXFwc0dHRLF68mKqqKgD+8Ic/MHz4cGJiYnjkkUf46aefWLduHb/73e+IjY3l+PHjHD9+nFmzZjFmzBiuuuoq0tK0YgSJiYn89re/ZerUqfz+979n5cqVLFu2DIDMzEymT59OTEwM06dPJysry+I5DVm5ciXXX389c+fOJTIykpdffpnnn3+euLg4xo8fz/nz2gizKX+++OILrrjiCuLi4pgxYwZnz2raKU888QSLFy8mISGBqKgoli9fbuf/Cdfk2+Pf6hQMhRDcMOwGJ3lkH/x8/EiISNDZkzKTLFZzUrgWLikTAPDkFwc5lFNi02sO79uZx+eOaLZNeno6b775JldeeSWLFy/m1VdfZdmyZSQmJrJp0yYGDx7MHXfcwWuvvcYdd9zBp59+SlpaGkIIioqKCA4OZt68ecyZM4cFCxYAMH36dF5//XUGDRrEzp07uf/++/n+e23Uc0lG18fHh5UrV9b7cUl698477+Stt97iV7/6FZ999pnunMYcOHCAlJQUKisrGThwIH//+99JSUnhN7/5De+++y4PPfQQS5YssejPpEmT2LFjB0II/vOf//CPf/yjXt43LS2NzZs3U1paypAhQ1i6dClGo2sKVQFEdY2y6fVKqkrYeGKjzj4hdAJ9g/ra9F6uQEJEAt8e/5aaupp6W2VtJduytjEjaka97eObPnaGe4pmcNmg7iway+4uX76cmTNn6qR3X3nlFZYtW1YvvTt79mzmzNHvymsovXuJS6N8aF5695NPPgE06d3//d//bfEcgKlTpxIUFERQUBBdunRh7ty5AERHR7Nv375m/cnOzubmm28mNzeX6upqIiMvb96ZPXs2/v7++Pv707NnT86ePUtoaGgLvek8bK2MuP7YerMAB9r889whc216H1chyD+ICaETSMpMMrNvPLGRqRFT8TFonz93zMn3dFoM6kKIMOBdoDdgAlZIKV9s1CYB+Bw4edH0iZTyqfY41tKI2l40XuwSQtCUPs4l6d1NmzaxevVqXn755foR+CUaSu9awtbSu/7+/vXvDQZD/e8Gg4Ha2tpm/XnwwQf57W9/y7x589iyZQtPPPGExev6+Pi4/Hx+fnk+YJugU1hRyNaMrTr71MipdOvQrd3Xd1VmRM3gh6wfzD7/hRWFpJ5JZUzfMQCsTF0JQGJsohM8VFjCmjn1WuBhKeUwYDzwgBBiuIV2P0gpYy++2hXQnUlj2d1JkyYxdOhQj5Hebc6f4uJi+vXrB8A777xjk/s5i4LyAgrKC2xyrfXH1uvkAPx9/blmwDU2ub6r0qtTL4tFqrdkbKl/vzJ1ZX1gV7gGLQZ1KWWulHLvxfelwGGgn70dcxaXZHdjYmI4f/48S5cuJSAggLfffpuFCxcSHR2NwWDgvvvuo7S0lDlz5hATE8OUKVPMpHefffZZ4uLiOH78OKtWreLNN99k1KhRjBgxgs8/1+c5N2b58uW8/fbbxMTE8N577/Hiiy+2eI61NOXPE088wcKFC7nqqqsICVGP1aB9OfyQ9YPOPi1yGkH+QU7wyLFMi5ymsx0pOEJOaY4TvFFYQ6ukd4UQEUASMFJKWdLAngCsBbKBHOARKeVBC+cvAZYAhIeHj8nMNNeUcLb0rpLd9RxsJb37/r73+SHTPKgH+Abwl+l/oaOfbStQOUt6tzmklDy59UlyS3PN7AkRCdwSfYuS3nUCNpPeFUJ0QgvcDzUM6BfZC/SXUo4CXgI+s3QNKeUKKWW8lDK+R48e1t5aoXAK+eX5/Jj1o84+I2qGzQO6qyKEYEr/KTr79uztKr3RRbEqqAshjGgBfZWU8pPGx6WUJVLKsovvvwaMQgi3e35XsruKhqw/th6TNJnZAo2BTI+a7iSPnMOEsAn4+/qb2apqq9iRvcNJHimaw5rsFwG8CRyWUj7fRJvewFkppRRCjEP7srDNKpVC0QYGdhvYrvMLKwr56dRPOvuMqBkEGgPbdW13I8A3gPGh43UZQFsytvDVrV95hDyCJ2HNSP1K4BfANCFE6sXXtUKI+4QQ911sswA4IIT4GVgOLJLOqpOnUAA+Bp/6XOq28O3xb6kz1ZnZOhg7WFw49AYs7TDNLc3ldOlpr/uSc3VaHKlLKbcBzX4VSylfBl62lVMKRXs5d+EcAD079mz1uSVVJWzL2qazT4ucRgdjh3b75o70DerL4O6DOVJwxMz+1NanmBg2kfvH3u8kzxSNUdovCo+ksKKwzaqCG45v0O0e9ff1Z3qkd82lN2ZKhH7BdHPGZlYfWO0EbxRNoYK6HVDSu3pmzZpFcHCwRSmFtpKYmMjHH9tWe6Ssuszi7tGEiASvyXhpitjesbo+kFJSUmVbjSZF+1BB3Q4o6V09v/vd73jvPdcvUPz9ye+prqs2sxl9jMyMmukkj1wHX4OvxUIgxZXFTUppKByPCuoNUNK79pPenT59OkFBze/AXL58eX1fLlq0SHdcSsmyZcsYPnw4s2fP5ty5c632ozkqayvZfHKzzj65/2Sv2D1qDZPC9XIVVXVVZBVnOcEbhSVcWqXx0m61htw04ibuH3s/5TXlXLtKX08xMTaRxNhE8svzWfDhArNj1ux6U9K7zpPe/dvf/sbJkyfx9/enqKhId/zTTz8lPT2d/fv3c/bsWYYPH87ixYttdv8fMn+gvKbczOZj8FGj9Ab0DepLVNcoThSeMLNvy9pG/+D+TvJK0RCXDurOQEnvOk96NyYmhttuu43rr7+e66+/Xnc8KSmJW265BR8fH/r27cu0aU2nF1qSB2iOWlOtRb308aHj6dqha6uu5elMCp9UH9TnDtY+X7tO72LhiIX4+fg50zUFLh7UmxtZBxoDmz0eEhjSJj0KJb3bNundTz/9lCeffBKA//znP8THNylN0SRfffUVSUlJrFu3jqeffpqDBw/i62v+EbXXRpcd2TsoqjR/OhBCeLwSY1uI7xvPmoNrqKq9PDiprK1kb+5ejyi+7e6oOfVGKOndtknv3nDDDaSmppKamtqmgG4ymTh16hRTp07lH//4B0VFRZSVlZm1mTx5MqtXr6auro7c3Fw2b9bPf1/iTNkZzpSdse7e0sS3x77V2eN6x9GrU6/W/UO8AH9ff8b2HQvAz2d/5uez2udn+6ntznRLcREV1BuhpHftI7171VVXsXDhQjZt2kRoaCjffmseROvq6rj99tuJjo4mLi6O3/zmNwQHB5OcnMzdd98NaF8cgwYNIjo6mqVLlzJlyuW86ccee4x169bV/15cWUxxZbFVvqWeSa3frNSQWQNnteWf6hVMDJsIQFZxVv0iaXpBepv3BihsR6ukd21JfHy8TE5ONrMp6V2FrbBWeldKyd+2/Y2MogyzNsN6DOOh8Q/Z3c+GuKL0blNIKXls82O8mfImcHlu/fqh1/M/g/7Hma55PDaT3lUoPJGj54/qAjqoUXpLCCEszp9vz96uctadjArqDVDSu97HhuMbdLb+wf0Z0r112TPeiKWgfrbsrMUvSYXjUEFd4ZEIIVrMlMkpzWH/2f06+9UDrlZyslbQPbA7nf0742swz1BSOuvORQV1hUcyuPtgBncf3Gyb745/p7OFBIYwus9oe7nlcbx13Vv8z0DzOfRdp3fpCnUrHIcK6gqvpKiyiJ2nd+rsM6JmYBDqz8JaRvcZjdHHfGdxeU05+87uc5JHCvXpVXgkOaU5zVa8//7k97oiGB39Otan6ims49kfn+V0yWmdfdfpXU7wRgEqqNsFJb3bMmlpaUyYMAF/f3/++c9/2uy6ERER5OfnU1pVSmlVqcU2UkqSMpN09qkRU3W1OBXNs+nkJrJLs3X2/Wf363R0FI5BBXU74O3SuxkZGSQkJDTbplu3bixfvpxHHnnEMU41oKquioqaCjObr8HXYsk2Rct0NHbUqVjWmmpJyU1xkkfejQrqDVDSu/aT3m1Mz549GTt2bLNKj3V1dSQmJjJy5Eiio6Prd+w2pKCggKuvvpq4uDjuvffeFnOkpZRU1Fbo7BPCJih53XZwSTagIWoKxjm4tKCXkt71LundxqSmpnL69On6vQOW5HiffPJJJk2axGOPPcZXX33FihUrmr1meU05JpNJZ58RNcM2Tnsp4/qN4/uT5mJ26QXpFFUWERwQ7CSvvBOXDurOQEnvtl1694YbbuDkyZNUV1eTlZVFbGwsAL/+9a/55S9/adHf5oiKiuLEiRM8+OCDzJ49m6uvvlrXJikpqb6fZs+eTdeumkyuj0HfP02VXovpFUPvTr1b7Z9Cy1UHiAiOoEfHHuRdyKs/JqUkOSdZfWE6GJcO6kp617Jfrii9C5r8LmjTWImJiWzZssWqf1tTdO3alZ9//plvv/2WV155hQ8//JC33npL187SRqGB3QbqbFV1VbpSdQAzB6giGG1l7U1r69+P6zeOr458ZXZ8Z/ZOFdQdjJpTb4SS3m2b9K49yM/Px2QyceONN/L000+zd+9eXZvJkyezatUqAL755hsKC5tWCbQ0Su8f3J9B3QbZzmkv5op+V+hsWcVZVksgK2xDi0FdCBEmhNgshDgshDgohPi1hTZCCLFcCHFMCLFPCOG2W/KU9K59pHcbc+bMGUJDQ3n++ed55plnCA0NpaTEPOiePn2ahIQEYmNjSUxM5K9//SsAr7/+Oq+//joAjz/+OElJSYwePZoNGzYQHh4OQHZJNtOunkZOjparXlNXo8t4AZgZNVNJArSDP278I3/c+EcAenXqZbGknVowdSwtSu8KIfoAfaSUe4UQQcAe4Hop5aEGba4FHgSuBa4AXpRS6r+2G6CkdxX2pLH0bkF5AWXVWtGNU8dP8V7Oe3Tr0I3/m/5/LrOD1J2kdy9xKZnh0lTnxhMb+ejgR2ZtenXqxZMJT6ovTxvRbuldKWWulHLvxfelwGGgX6Nm1wHvSo0dQPDFLwOFwunUmeq4UHNBZ58WOc1lArqnEN83Xhe8z5adJbtEv0FJYR9a9YkWQkQAcUBj0Yx+wKkGv2ejD/wIIZYIIZKFEMl5eXmNDzsdJb3rmZRVl+kWuwN8A5gUbpt1CsVlggOCLS5S787Z7QRvvBOrg7oQohOwFnhIStl4xcnSc5VuXkdKuUJKGS+ljO/Ro4fF+yiBfYUtaZzGKKVEIpkUPokOxg5O9Mxzie+rnxlIzklWf9sOwqqgLoQwogX0VVLKTyw0yQbCGvweCjStptQEAQEBFBQUqP98Rbsx+hgx+hi5UHMBk9Q2G0kpKS8p53zNeaZFTnOyh55BaOdQQjub71cY02eMbgqmoLxAFc9wEC3mqQvtf+dN4LCU8vkmmq0DlgkhVqMtlBZLKXNb60xoaCjZ2dm44tSMwj3ZX7m/Xo1RIjlfc54CY0H9phlF+3h//vs6W5B/EENDhnI477CZfXfObiK7RuraK2yLNZuPrgR+AewXQlzasfIoEA4gpXwd+Bot8+UYUA60fvsgYDQazXYxKhTtIS0/jZXbV+rsv5/0e31jhU0Z23esLqjvydnDwuELVRaMnWkxqEspt2F5zrxhGwk8YCunFIr28tD6h/j5zM+66kdRXaOI6hrlJK88j4fWPwTAC7NeMLPH9o5l1f5VZpr1RZVFHDt/jEHd1WYve6LyuRQeya7Tu0gvSNfZlSSAbUk9k0rqGb3kREe/jgzvMVxnV1kw9kcFdYVHUliplwvoHtid2N6xTvDGO7Ekx7s3d2/9wrXCPqigrvA4yqrLLOq8TI+crjYbOZBRvUfhazCf4S2tKuVowVEneeQdqE+4wuNIykyyuNnoyvArneSRdxLgG8CIniN09j25e5zgjfeggrrCo6g11bL55Ga6+Hehi3+Xevuk8EkE+AY40TPPZHD3wbrF6IZY2oikpmDsi0vrqSsUrWX36d2UVJUwuf/kepsQQm02shMr5jZfaSqmVwy+Bl9qTZf190urSjlScIShIUPt7Z5XokbqCo9BSsnGExt19jF9xqjNRk4iwDeA6F7ROvueHDUFYy9UUFd4DOkF6fVqgEmZSSRlJgGq/qg9WfLFEpZ8saTZNmP6jNHZ1BSM/VDTLwqP4bvj3+rUIfsAACAASURBVNW/L64qBmBAtwFqa7odOVJwpMU20b2iMfoYqamrqbeVVZepKRg7oUbqCo8gtzSXA+f0sskzo9RmI2cT4BvAyJ4jdfbknGQLrRXtRQV1hUew6eQmnc1oMDKq9ygneKNojKUpmJTcFDUFYwdUUFe4PaVVpWw/tV1n79qhq9ps5CLE9IrB6GM0s12aglHYFvWJV7g9WzK2mKXMgVYXU1U2sj+xvWOtkl7w9/W3OAWjsmBsj1ooVbg1NXU1bMnYorP/efKfmT9svuMd8jIaqzM2x5g+Y0jJTTGzpZxJ4ZboW9QTlQ1RPalwa3Zk76CsuszMZhAGpkZMdZJHiqaI7hWttGAcgArqCrdFSsl3J77T2cf2G8uD3zzI7Z/c7gSvvIvbP7nd6n5uSgtmb+5eW7vl1aigrnBb9p/bz9myszr7zKiZZJdk129EUtiP1vaz2ohkf1RQV7gtDTcbXWJoyFDCuoRZaK1wBSzJ8ZZUlXD8/HEneeR5qKCucEsyijIspsOpykaujZLjtT8qqCvcEkuj9D5BfRjRQx8wFK7F6D6jdba9uXt1GviKtqFSGhVuR355vsWR3cyomfWV6ieETnC0W15JW/p5VK9ROjne4spijhceZ2C3gbZ0zytRQV3hdmw8sVE3qusS0IUrQq+o//2vM/7qaLe8krb0cwdjB4b3GM6+s/vM7Hty9qigbgPU9IvCrbhQfYEfs37U2adFTtMtwClcFzUFYz9a/CsQQrwFzAHOSSl1+3yFEAnA58DJi6ZPpJRP2dJJhXdSWwuVlVBTA9XV2s+NmVspKKqm4d++n48/YabJZGaCjw8YjXD3dzdiEPDxTWvx9YWLszIKG3PjhzcCsPamta06b1TvUfgYfKgz1dXbiiqLOFF4ggHdBtjUR2/DmqHNSuBl4N1m2vwgpZxjE48UHkt1NRQWaq+ePaFbt6bbfvoprF9vbjNRw67O31PdKED3q5rE8s2BZrY93QsAWLYMHn8c+vZt+l5paWAyQXAwdO0KHTq05l/l3RSUF7TpvEBjIMNChunkkvfm7lVBvZ20GNSllElCiAj7u6Jwd6SEsjI4cwbOnYO8PO1nfj4UFGjHLrFoEUxtZie/v7/edsbvJ6pFqZlNYKBf1fRm/Qpood70J59AZubl3zt0gO7dtVfPntqrRw/o3VsL/GrUbxvG9B2jC+p7cvewYPiC+gVvReux1STkBCHEz0AO8IiU8qClRkKIJcASgPDwcBvdWuEMKivh1Ck4ffry68wZuHDBuvMLC5s/3jgQS0xk+2/QtetRPZYA2Xz90ZaCelGR+e8VFZCdrb0a4++vBfd+/S6/wsKgU6fm76HQM6rXKAzCYLabtLCikIyiDFWtqh3YIqjvBfpLKcuEENcCnwGDLDWUUq4AVgDEx8erFRE3o7QUVq+GrCxtBN4eGgfSxjQOxPnGvVQa8nXtQquubvFezQV1kwlKSlq8RD1VVdqovuHIHrRpm3HjYL4ShrSajn4dGRoylEN5h8zse3L3qKDeDtod1KWUJQ3efy2EeFUIESKl1P8FKtyaDh0gNVVbwGwv1gR1g0H7afSTpPuvp6OPNvVx6dXbMJJRIaFIqQXnujrtVVsLA0zTqTNBYKB2naYoLgZbJFwUFmojfG9jemTzU18tMabvGH1Qz9nDjcNuVFMwbaTdQV0I0Rs4K6WUQohxaGmSbVs9UTiFoiJtsfDoUbjlFvBt4lPh6wvh4XDiRNvuYzBoc9LBwRAa2nzbuDh49VUteB/OS6N4xyldm4cnXsPgJmde/my1X2PHan1waRG3rq7lcyzR0ozizz/Drl0wZAgMHarN07t73PrzFOv72RKxvWNZtW+V2RTM+YrzZBZnEhEc0U7vvBNrUhr/CyQAIUKIbOBxwAggpXwdWAAsFULUAhXAIqmSTV2a6mo4cgQOHoTDhyE39/KxK6+EqKimz42MbD6o+/pCr16XXz16aK/u3bUpiuZGzQ1pGOzWH1uvOx7ZNZJB3SzO8rWKrl3h7rsv/y6lNnovKNAWevPy4OxZ7XXmjNZ3TdFSUN+/H5KTtRdo2T/DhsHIkVqQ90Y6+XViSMgQDucdNrPvydmjgnobsSb75ZYWjr+MlvKocGEKCrSgsm8fpKc3PYVy5EjzQT0i4vL7nj21RcJLC4Z9+0JIiPWB2xoyijJIy0/T2a8ZcE2zj+f/s+p/APjmtm9adT8hLj9NDGiUWSelNpLPybm8OJydfflLsV+/5q99pJH+2Pnz8OOP2stggFPdtQXXnBzo08c9RvFt7eeGjOkzRh/Uc/cwf9h8NQXTBtQWPA9FSi3gpKZqL0uZHJY4cgRmzWr6+PDh8PDD2qi0pawSW/DNUX2w6NWpF6N6j2r2vIoa209wC6GNrrt100bXl6ip0UbyTU1bgTb6P6uXfq/HZNLm5Csq4MkntS/H2FhtGioqyrZflLbEFv0c1yeOD/Z/YDYFU1BeoKZg2ogK6h6ElFpmyp492iu/DUvVx45pAaapINKpEwwe3D4/rSWnNIfUM6k6+6yBs1yqpqXR2PIaQeNRekvk58PGjdorKAhGj4YxY2DQINcN8G1FTcHYFhXUPQAp4csvYceOtgVy0KZRhg51rbldS3PpXTt0ZVy/cU7wpn3ExmpPOOnp2jrGyZPal6c1lJbC1q3aq3NnuP56be3Dk1BTMLZDBXUPQAhthN2agB4UBCNGXF6kCwqyn39tIb88n92nd+vs1wy4xi2Fu4xG7Qln8GCYO1fbvJWeri1WHzigrXlYQ0mJ5d227k5s71g1BWMj3O+vQ2GRK67Q0hKbIywMRo2CmBhtTtyVB0Abjm/Q1a0M8g/iynDrhqhzBru2FFFAgPZ/MWqU9qR19izcsVKTUhAXms6d9/fX/v9cBVv1c5B/kJqCsREqqLsBl7bfN87GaEhcHKxaZZ7VIoQ2BxsXpz3+Nyeg5UoUVRZZlNedHjkdPx8/q67xyMRHbO2W3RBCkx7o2lV7/fMPWpZSSgocOmT+fxoXB37NdIGU8MUX2hNYZKT9v7ht2c9qCsY2qKDuotTWalkrW7dqi2z9+8OjjzbdvkMHiI7WzhkyRFtUi4tzvWkVa9hwfINZVRzQalsmRCQ4xyEH06kTTJyovSortQCfnKxN1YxrYTkhKwu++kp7hYbClCnaOY7IVGovagrGNqig7mIUFUFSEvzwg7kmSWamJqAVFtb0ufPnaztCu3Sxv5/2oqSqhKTMJJ19auRUOhit18RNWJkAwJbELTbyzDkEBGhBedw4KC9vOThv23b5fXa29vS2di1MmKCpYvbqZVv/bNnPTU3B7D69WwX1VqCCuotw4gRs2gR79zadFfHjj5pkbVP07Gkf3xzJd8e/o6auxszm7+vPjKgZTvLIdQgMbP54dbUmQ9CYykrYvFl7DR8O06dri+SuOKMR3zfe4hSMkuO1HhXUnYjJpM2bbtxonZ7Kzp1w441aJoUnUlpVypaMLTp7QkQCnfyUtm1LpKRoAbw5Dh3SXr17a8F9wgTX+jzF9Y7TacEUVhSqikitQAV1J1BdDT/9BN99Z30aYp8+cNVVtlEUdFU2nthIdZ25uIrRx8jMqJlO8si9iIuDxERt+q6lQcKZM9rUzLp12rRMQgJ07OgIL5uno19HhvUYxsFz5iUZdufsVkHdSlRQdyDl5bBlizbN0rAKUFMYDNpOwoQEGDjQNR+XbcWF6gtsztiss0/pP4Ugfzdc7XUCfn7ayHvCBG0+fetWbUNacyJkpaVaYF+/HiZPhpkzNd0bZxLfN14X1Pfm7uWmETe51E5iV0UFdQdQWqpNsWzZ0vLjMWi7BidP1l7uvOjZGr49/i1VtVVmNl+DL1cPaLkIhiVuGnGTLdxyW0JD4bbb4IYbYPt27bPXXGGT6mrtM7p5s/alMHeudcHdHv0c2zuW9w3vmxWlLq4s5tj5Ywzu7iCNCjdGBXUHkJSkL6JsibAwmDED4uObF4fyNEqrStl8Uj9Kv6r/VXQJaNu32v1j72+vWx5BYKA2dz5tmqbSuXGjtpO1KerqtNH93LnWXd8e/RxoDGREjxHsO7vPzJ6ck6yCuhV4UehwHlOnwoYNTY/SR46Eq6/WtpB78hRLU3x7/FuLc+mzBjYjF9kC5TXlgBYgFNrnKiZGe506pa3n7N5tOdNq0iTrp2Ds1c/xfeN1QX1v7l4WjVykpmBaQAV1B3BptPTVV5dtBoOWe3zNNZoOubdSXFlsMeNlSv8pBAe0fXL32lXXAu6fp24PwsJg8WJNGGzjRm1PxKV5d1/f5qWXG2Ovfh7VexS+Bl+zTWilVaWk5acxvMdwm97L01BfeTagulorSNwc06drG0d8fbWFz2eegV/+0rsDOmij9MZ56UYfI9cMvMZJHnkP3brBTTfBX/8Kc+Zog4+JEzWpguY4fty6hf72EOAbQHSvaJ1912kLifgKM9RIvR2YTNqGoC++0B5Z581rum3HjnDPPdooyVsWP1uisKKQrRlbdfapEVPp7N/ZCR55J506aXPoM2e2XFS8qgpef10byFxzTfPa++1lXL9xpOSmmNlSclO4Lfo2jD4ulFzvYqig3gak1PQ4PvlEy/cF7TE2IUHLXGmKhtVyFPDV0a90Gi/+vv5tznhRtA9r9GG+//6yfMXnn0NGL63+rD2Ce3TPaAJ8A6isvbwYVVlbyYFzB4jrE2fbm3kQavqllWRlwfPPa5XuLwV00EYw37S9TKPXcbbsrEUlxmmR01ReuotSXq4t+DektlaTDX76aU1wzJYYfYzE9o7V2dUUTPOokbqVFBXBZ59p6V5N7epMStJSErt3d6xv7si69HU6vfRAY6DNRumJsYk2uY7iMhs2aIG9IYPLEwHIKYDly7Wn0QULtB3QtmBcv3HsyN5hZtt/bj+VtZUE+LqB9KQTUEG9BWpqtKmVb75peTF01CjvTElsLVnFWSTnJOvs1wy8xmapcSqo255p07TC2ElJl1Mhh1QkmrU5cEDTlpk69fLia3sYGjKUIP8gSqtK6201dTWknkllfOj49l3cQ1FBvRn27YM1a1rWZxk0SBudREQ4xC2357O0z3S2zv6dmRY5zWb3yC/X/tNCAkNsdk1vp3NnTdp5+nTtqXXPHqg0aP0cYLrczyaTJoWxc6eWNjlpUtsHOz4GH8b0GaNLe911epcK6k3QYlAXQrwFzAHOSSl1S31C08N8EbgWKAcSpZR7be2oI8nL04L5/v3Nt+vVSwvm0dFqhG4taflpOl0PgNmDZ1td1cgaFny4AFB56vagZ09YskQTDZvx/gIqK2FuwRZdu7IyeP99TeP9llvaPugZ22+sLqgfzjtMSVWJypKygDULpSuB5rYj/A8w6OJrCfBa+91yDjU1WnriE080H9A7dtR0zR9/XNuhpwK6dUgp+fjQxzp7SGAIk8InOcEjRXuIitJq3fbp0/w6UkYG/O1vWoC/cKH19xnQdQDdOpjXYjRJk8XC5AorgrqUMgk430yT64B3pcYOIFgIYaNlEsdx6BA89RR8+WXTuboGgzZX+PTT2k8fH8f66O7sOr2LU8WndPZ5Q+bha1Azge5KUBA8+SRcd13T9VOl1HauPv64VgimNQghGNdPX8dv5+mdbfDW87FFSmM/oOFfavZFmw4hxBIhRLIQIjkvL88Gt7YNUmryo82p2A0dCn/+szZCdwXdaXejpq7G4lx6eJdwi3+wCvfCaIRrr9UGPFdc0XS70tK21QSwNH+eWZTJmbIzFlp7N7YI6pYmHyz+t0kpV0gp46WU8T169LDBrW2DEJpMqaVplK5d4d574aGH1Jb+9vD9ye85X6F/4FNlyjyL4GBNV+aRR6CfhaHdiBFajYDW0ieoD+FdwnX2xumOCttkv2QDDcshhwI5NriuQwkL01b1N27UfjcYtG3Ts2eDv79zfXN3yqrL+Pro1zp7TK8YhoQMscs9l8Yvtct1FeY01c+DBsGf/qTps69bpymUGo3agmlbv8PHh44nqzjLzLYzeyfXDblODQwaYIugvg5YJoRYDVwBFEspc21wXYczdy4kJ0OPHtrI3VYbKLydz9M+N9vqDdo86fxh8+12z5tH3my3aysu01w/GwzaQGnMGPjwQ61wR3se0Mf2G8vHhz4227R2vuI8R88fVTrrDbAmpfG/QAIQIoTIBh4HjABSyteBr9HSGY+hpTT+0l7OtoczZ7QCAJYeCS8REAC//7025aK++G1Ddkk2P2T9oLNPCp9EnyD7fWteWpAN6xLWQktFe7Cmn4ODtRRIa+bSP/pIq7U6cKD+WGf/zgzvMZwD5w6Y2Xdk71BBvQEtBnUp5S0tHJfAAzbzyMaYTNqUyuefa3nljz7afFWhbt2aPqZoHVJK1hxYg2z01xzgG8C8Ic1IWtqAX3z6C0Dlqdub1vRzSwOlvXu1v9VNm7QRvqVsmvGh43VBfU/OHm4ZeYtSbryIRwt6nTsH//wnrF2rpSmePm1eqEJhX/bm7uVIwRGdffbg2WrTiMKMCxfggw+091Jqwf3pp7UNTg0Z1XuUTvOlsraSlDPmEr3ejEcGdSm1SupPP60J+jdk/XrIzHSOX95ETV2NxY1GvTr1sqkcgMIz+PBDLd2xIefOwT/+oT1lX9o74ufjx+g++vQZS4qf3orHBfWiInjpJe1bv7paf9xkUqN1R/D10a8tpjDeNOImtdFIYUZNTdN7RKSEr7/WdqTmXMypuzL8Sl27tPw0CsoL7Oil++BRQT0lRdsV2pSusxDaXN3ddzvWL2/jbNlZNhzfoLOP7DmSkT1VpRCFOUYj/O53cOONTa93nToF//d/WpGOqOAB9OzYU9fmp1M/2dlT98AjhkyVldrj24/NPIGFhEBiopY/q7AfUko+2P+BrqKRj8GHm0bc5DA/Hp7wsMPu5c3Yqp8NBrj6ak1L6e23Nb2YxtTWakJ7Bw4IYqdcyYYLn5od/+nUT8wZPMfrc9bdPqhnZsJ//tP8Fv/JkzU1RbWJyP4k5ySTlp+ms18z4Bp6derlMD/mDpnrsHt5M7bu5969tbTib77RdJhMJn2bgwchPXM8haM+o2vXy5lV5yvOk5afxrAew2zqk7vhtkFdSvjuO03Xua7OcpvOneHOO1VtUEdRUVPBhwc/1NlDAkO4dtC1DvUlPT8dwG47VhUa9uhng0HbyT1yJLz1lnnZyEvUlgVzOnUE5/sfIDLicn3UH0/9qIK6sx1oC6Wl2iNaczURY2PhF7/QKqUrHMMnhz+hpKpEZ7955M0OzyG+98t7AZWnbm/s2c/9+8P/+39agffNm/XHe1dfyaHTByguhmFDoUMHSMlN4UL1BTr6ea/qntstlKanN1/k1s9PC+b33acCuiM5UnCEpMwknT22dywxvWKc4JHCE/Dz05RRH3xQe/JuSLeaGIyyE2VlsDcF8vKh1lTr9SJfbhXUpdQ2EhUXWz4eHq6JCLWnfJai9dTU1fDez+/p7P6+/kqDRWETRo6Exx4zn0o14Euvak2St66Oem3YrZlbdbuYvQm3CupCwF13WV7wnDlTW2Dp5bi1OMVFvjjyBecu6Feqbxh6g65ijULRVoKCYNkyuOmmy6mPvasnA9Cn92WxsLNlZy3uZPYW3Cqogxa0b7vt8u+dOmmPZgsWNK/porAPmUWZFnPSB3QbQEJEguMdUng0l/aa/P73Wq3UQFMvwjsOIWqAebutmVud46AL4JZh8IortLn1vDxt5B4c7GyPvJOauhreSnlL96jra/DljlF3ODVf+E+T/+S0e3sTzurn8HBtEXXNGrghdgqfZqWbHU/JTfHawtRuGdRBE9v38bmcyqRwPJ+mfWqxnNicwXPo3am3Ezy6zIyoGU69v7fgzH4OCNBSlmtNo9h0trNZ5pVJmtiWta0+lbaiQsuO8QbcNiQajSqgO5O0/DQ2ndiks4d1CePqAVc7wSNzUs+kknom1dlueDyu0M++Bl8mhU/S2X/I/AGTNLFjh7bIevSoE5xzAiosKlpNeU05K1NX6uy+Bl/uirsLH4OP451qxEPrH+Kh9Q852w2Px1X6eVL4JN103/mK82zct5/334eSEnj+eS3f3dMTY1RQV7QKKSXv73ufwopC3bEbht1g12pGCkVTdA/srhOLq62FZ9dupKZG+91kgtWr4Z13qLd5IiqoK1rFtqxt7MnZo7MP7j6Y6ZHTneCRQqExNWKq2e8nT0J25RHKDKfM7Nu3a8Vziooc6Z3jUEFdYTWnS06z5uAanT3AN4DE2ESvV8dTOJfhPYabPSn276/tQj3tv1HXNiNDk/JtXETHE1BBXWEVVbVV/Hvvv6mp0z+33h5zO90DuzvBK4XiMkIIs6dFPz9NytcYsZtqod+Gfmme/ScPk2F325RGheOQUrJq/ypyS3N1xyaFT2Jsv7FO8Kp5/jL9L852wStwtX4eHzqez9I+o6y6DACDgKgBdURFbCHru+vqy+JdorZWm2PPyYH58z0jo84D/gkKe7M5YzM7s3fq7H2D+rqstsvEsIlMDJvobDc8HlfrZ6OPkcn9J+vsZzts5cGHqunSxfJ5330Hr7yiFdxxd1RQVzTL0YKjfHTwI53d6GPknjH34Ofj5wSvWuanUz+p8mYOwBX7OSEiQZdWe6H6AmeMP/HooxAZafm8Awfg73+HAjcvdWpVUBdCzBJCpAshjgkh/mDheKIQIk8IkXrxpaqAegCFFYW8secNTFJffua26NvoG9TXCV5Zx6ObHuXRTY862w2PxxX7uUtAF8b21U8Jrj+2nk6da3n4YU1qxBI5OfDXv8KJE3Z20o60GNSFED7AK8D/AMOBW4QQwy00XSOljL34+o+N/VQ4mKraKl7Z/QqlVaW6YwkRCUwIm+AErxQK67C0q7mwopCd2TsxGuGXv9Tm0C0lbJWWwnPPwR595q5bYM1IfRxwTEp5QkpZDawGrrOvWwpnYpIm/rP3P5wqPqU7NqDbABaOWOgErxQK6+nXuR+xvWN19m+OfYNJmhACrrkGli61LOVdWwsrVkCqGypNWBPU+wEN/7qzL9oac6MQYp8Q4mMhRJilCwkhlgghkoUQyXl5eW1wV+EIPj70MfvO7tPZgwOCuS/+PnwNKmlK4fpYqoubdyGP3ad31/8+ahT87neWlV4jImC4pTkJF8eaoG5pR0lj9YQvgAgpZQywEXjH0oWklCuklPFSyvgelxTtFS7FphObLAp1GX2MLB271CulTBXuSf/g/ozoOUJn/+bYN2Zy0WFh8Mc/apuVLtG9OzzwgJbr7m5YM+TKBhqOvEOBnIYNpJQN14v/Dfy9/a4pHM2O7B18ePBDnV0IwV1xdxERHOF4p9rIC7NecLYLXoGr9/O1g67l4Dnzgsa5pbkk5ySb7a8IDoaHH4Y334QjRyzXRHUXrAnqu4FBQohI4DSwCLi1YQMhRB8p5aWdKfOAwzb1UmF3fj7zM++kWnzAYv6w+cT1iXOwR+3D0nyqwva4ej8P7DaQwd0H68rbfZ7+OXF94symEv39tYL1eXnuXRazxekXKWUtsAz4Fi1YfyilPCiEeEoIMe9is18JIQ4KIX4GfgUk2sthhe1Jz09nxZ4VFlMXr+p/FTOjZjrBq/ax8cRGNp7Qa34obIs79PPswbN1trwLeWzL2qazGwzWBXRXlu+1asVLSvk18HUj22MN3v8R+KNtXVM4grT8NF7e9TK1plrdsbg+cdwafatbCnU9k/QMoCog2Rt36OehIUMZ1mMYh/PMJxC+PPIl40PHE+Ab0KrrlZfD8uUwZw6MHNlye0ejdpR6MYfzDvPyrpctinQNDRnK3aPvxiDUR0Th/swfNl9nK60qbfVTRk0NvPqqJuv7yiuajK+rof5ivZSD5w7yyu5XLAb0yK6R3D/2fpW6qPAYwruEM67fOJ19w/ENZrVNm8NkgrfeulwWz2SClSthwwYbOmoDVFD3QnZm72xyhB7eJZwHxz2Iv6+FHRkKhRtz3dDrdJowVbVVrD201qrzv/sO9u7V29eu1V6uMs+ugrqXseH4Bt5KecviomhEcAS/mfAbOvp1dIJnCoV9CQkMYUr/KTr7juwduuwYS0yZAkOHWj62YQO89542enc2Kqh7CSZpYvWB1U2OSiK7RvLr8b8m0BjoYM/swxtz3uCNOW842w2Px936+dpB11r8jH+w/wOLyQINCQjQ8tfj4y0f//FH+Pe/0Wm2OxoV1L2AC9UXeHHHi2w+udni8cHdB/PrKzwnoAMMCRnCkJAhznbD43G3fg7yD7K4aJpbmst3x79r8XxfX7jrLm3Ubom9e7UF1Kqq9nradlRQ93BySnP467a/kpafZvH46D6j+fX4X9PB2MHBntmXL9K/4Iv0L5zthsfjjv08KXwSUV2jdPavjn5Ffnl+i+cbDHDLLVpKoyUOHYIXX4SKivZ62jZUUPdQpJT8mPUjf/nhL+RdsCyelhCRwD1j7vHILJfntj/Hc9ufc7YbHo879rMQgttibtOl69bU1fDm3jctrjfprwFz58JNN1k+fvw4/OtfUFZmC49bhwrqHkhlbSVvp77Nuz+/azHDRQjBwhELWTRykcpDV3gloZ1DmRY5TWc/UXiCr49+beEMy0yfDomJlnXZMzM1XfZifc1ru6L+oj2MtPw0ntr6lMWaogAdjB14cNyDzIia4ZY7RRUKWzF3yFy6B3bX2b888iXHzx+3+joTJsC992rz7Y3JyYF//hMKC9vjaetQQd1DqKytZNW+Vfxr+78oKLdcZLFPUB/+OOmPFuVIFQpvI8A3gMVxi3WDGyklb6a8SUWN9ZPicXGaVK/RqD927hw8+yzktzxdbxNUUHdzpJTsPr2bxzc/TlJmUpPtJoZN5I+T/kivTm4sP6dQ2JiB3QZaLKZRUF7Av/f+26r59UsMHw6//rWW+qi7XgE8/7wmM2BvPG+FzIvILslmzYE1zW6c8Pf157bo27giP7wqPwAACqZJREFUtIlKux7Keze852wXvAJP6Oc5g+dwOO8wJwrNq00fPHeQNQfWcEv0LVZfa9Ag+O1v4YUXNOEvs/vMsTyStzUqqLsh5y6cY136OrOyXJYY1H0Qd466kx4dva/KVFgXixUVFTbGE/rZIAwsjlvMM0nPUFlbaXZsS8YWenfqzdTIqVZfr39/reBGw+yX226DiRNt6XXTqKDuRpwtO8u3x79l+6ntzT4W+vn4MX/YfBIiErx2MXTNgTUA3DzyZid74tl4Sj/36NiDu0bfxau7XzUrdQew5uAaOvl1MquU1BKhofDII1pgv+YamDzZ1h43jQrqLo6UkhOFJ9h4YiMpZ1J0H7jGjOo9iptG3ERIYIiDPHRNXkt+DXD/YOPqeFI/x/SKYeHwhbqSjpcWTmtMNUwMs3643acPPPEEBDp4o7YK6i5KZW0lu07vYmvGVrJLslts37NjT24eeTMje7qgar9C4SZMi5zGmbIzuqQDKSXvpL5DdV01CREJVl/P0QEdVFB3KWpNtRzOO8zO0ztJPZNqceNQY4L8g5g9aDZX9b/KI3eGKhSORAjBopGLKKgo0BWsBvjv/v9y7sI55g+b77J/b67plRdRUVPBobxDpJ5JZf+5/VbnxgYaA5kRNYMZUTOU9rlCYUN8DD4sjV/KG3veYP/Z/brjm05s4kThCe4ZfY/FzUvORgV1B1NVW0VGUQZHCo5wKO8QGUUZrcqFDQ4IZuaAmUwKn9Tq2ooKhcI6jD5G7ou/jzf3vsneXH1ljJOFJ3km6RnmDZnH5P6TdcU3nIkK6nak1lRLbmkuWcVZZBVncbLoJKeKT7UqiF8iqmsUUyKmEN833mUf+1yJj2/62NkueAWe3M++Bl/uGXMP76S+w47sHbrj5TXlrD6wmi0ZW7hx+I1E94x2iWwzFR3aSa2plqLKIgrKCzh34Rx55Xmcu3COnNIc8i7ktSmAX6KjX0fi+8ZzVfhVHpEP7Ei8PfvHUXh6PxuEgcTYRMK6hLH20FqLf89nys7wyq5XCAkMYXL/yUwIm0Bn/85O8FZDBfWLSCmpNdVSVVdFVW0VFbUVVNRUUF5TzoWaC5RVl3Gh+gKl1aWUVJVQUlVCYUUhpdWlLaYZtoYA3wBG9hzJFaFXMLzHcDUqbyMrU1cCkBib6FQ/PB1v6GchBDOiZhDVNYoVe1ZQWGFZnSu/PJ9PDn/CJ4c/IaxLGMNChjGg2wD6dOpDj449HKaIalXEEELMAl4EfID/SCn/1ui4P/AuMAYoAG6WUmbY0tGKmgo+T/8cKSUSafGnSZp0rzpZp/001VFrqq1/1ZhqqDXVUl1XTU1dDdV11e0aVbeH7oHdGdFjBLG9YxkSMkQFchvgDcHGFfCmfo7qGsWfJ/+Zjw59xPZT25tte6r4FKeKT8FFsUcfgw/BAcEE+QXRya8TQf5BBPkFEdcnzmLBjvbQYvQQQvgArwAzgWxgtxBinZTyUINmdwGFUsqBQohFwN8Bm+5GqDHVNFmOzd3o1qEbA7oNYHD3wQwLGeaV2/gVCneko19HEmMTmRoxlQ8Pfsix88esOq/OVEdBeYFOQbVnx56OD+rAOOCYlPIEgBBiNXAd0DCoXwc8cfH9x8DLQgghbTgvIXD+AkRbCA4IJrRzKOFdwgnvEk5EcARdO3R1tlsKhaId9A/uzyMTH2H/uf1sydhiMafdGjr5dbKxZ9YF9X7AqQa/ZwONJf/q20gpa4UQxUB3wExBWAixBFgCEB4e3kaXXY/O/p3p1qEbIYEh9OzYkx4de9C7U2/6dOrjcbU/FQqFhhCCmF4xxPSKIb88n21Z20g9k0puaa7V1wjyD7K5X9YEdUtD5MYjcGvaIKVcAawAiI+Pb9Uo3hGLDAZhIMA3AH9ffwJ8Awg0BtLBtwMd/TrS0diRjn4d6eTXiS7+XegS0IXO/p0JDghWc+AKhZcTEhjC9UOv5/qh11NUWcThvMMcLzzOmbIz5JbmUlZtuVhpkJ9zgno20DCfLhTIaaJNthDCF+gCnLeJhxfx8/Hj5pE3IxAIIcx+GoQBIS7+bPC7j/DBx+Cj+2k0GPE1+GL0MeLn44fRoP10pQ0Eivbx9W3W15lUtB3Vz3qCA4KZEDaBCWET6m3lNeWUVpVSWl1KWXUZpVXaz+CAYJvf35qgvhsYJISIBE4Di4BbG7VZB9wJbAcWAN/bcj4dtB1elgrFKhSWCDQ6QUnJC1H9bB2BxkACjYH0wv6Vx1oM6hfnyJcB36KlNL4lpTwohHgKSJZSrgPeBN4TQhxDG6EvsqfTCkVLvLr7VQDuH3u/kz3xbFQ/ux5WTQZLKb8Gvm5ke6zB+0pgoW1dUyjaziVNbBVs7IvqZ9dDFZ5WKBQKD0IFdYVCofAgVFBXKBQKD0IFdYVCofAghI0zD62/sRB5QKZTbm6ZEBrtgFVYRPVTy6g+ahnVRy3TVB/1l1I2KRjltKDuagghkqWU8c72w9VR/dQyqo9aRvVRy7S1j9T0i0KhUHgQKqgrFAqFB6GC+mVWONsBN0H1U8uoPmoZ1Uct06Y+UnPqCoVC4UGokbpCoVB4ECqoKxQKhQfhdUFdCDFLCJEuhDgmhPiDheO/FUIcEkLsE0JsEkL0d4afzqSlPmrQboEQQgohvC41zZo+EkLcdPGzdFAI8YGjfXQ2VvythQshNgshUi7+vV3rDD+diRDiLSHEOSHEgSaOCyHE8ot9uE8IMbrFi0opveaFJh18HIgC/ICfgeGN2kwFAi++XwqscbbfrtZHF9sFAUnADiDe2X67Wh8Bg4AUoOvF33s6228X7KMVwNKL74cDGc722wn9NBkYDRxo4vi1wDdo1eXGAztbuqa3jdTri2hLKauBS0W065FSbpZSll/8dQdapSdvosU+usjTwD+ASkc65yJY00f3AK9IKQsBpJTnHOyjs7GmjyTQ+eL7Lugrqnk8Usokmq8Sdx3wrtTYAQQLIfo0d01vC+qWimj3a6b9XWjfkt5Ei30khIgDwqSUXzrSMRfCms/RYGCwEOJHIcQOIcQsh3nnGljTR08AtwshstHqNTzoGNfcitbGLOuKZHgQVhXIBhBC3A7EA1Ps6pHr0WwfCSEMwL+AREc55IJY8znyRZuCSUB72vtBCDFSSllkZ99cBWv66BZgpZTyOSHEBLTqaSOllCb7u+c2WB2zLuFtI3VrimgjhJgB/D9gnpSyykG+uQot9VEQMBLYIoTIQJvnW+dli6XWFmP/XEpZI6U8CaSjBXlvwZo+ugv4EEBKuR0IQBOxUlzGqpjVEG8L6vVFtIUQfmi1VNc1bHBxauENtIDubfOg0EIfSSmLpZQhUsoIKWUE2rrDPCllsnPcdQotfo6Az9AW3RFChKBNx5xwqJfOxZo+ygKmAwghhqEF9TyHeun6rAPuuJgFMx4ollLmNneCV02/SOuKaD8LdAI+EkIAZEkp5znNaQdjZR95NVb20bfA1UKIQ0Ad8DspZYHzvHYsVvbRw8C/hRC/QZtSSJQXUz68BSH+f3t2TIAwFANA9BhwgCpcoAQbdYCp6uj6GcrOnr4nIEOGG5Lbp/NE9/j9Ft7VvWqttXX+Gp7VXh3V6+/Mi+0QYLSrnV8ARhN1gEFEHWAQUQcYRNQBBhF1gEFEHWCQL5VwLSGFolwyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prior_alpha, prior_beta = 2., 2.\n",
    "posterior_alpha, posterior_beta = prior_alpha + 7, prior_beta + 10 - 7\n",
    "\n",
    "posterior_mean, posterior_var, posterior_skew, posterior_kurt = scipy.stats.beta.stats(posterior_alpha, posterior_beta, moments='mvsk')\n",
    "\n",
    "xs = np.linspace(scipy.stats.beta.ppf(0.00001, posterior_alpha, posterior_beta), \n",
    "                 scipy.stats.beta.ppf(0.99999, posterior_alpha, posterior_beta), 100)\n",
    "plt.plot(xs, scipy.stats.beta.pdf(xs, prior_alpha, prior_beta), 'b--', lw=5, alpha=.6, label='prior')\n",
    "plt.plot(xs, scipy.stats.beta.pdf(xs, posterior_alpha, posterior_beta), 'g-', lw=5, alpha=.6, label='posterior')\n",
    "plt.axvline(posterior_mean, label='posterior mean')\n",
    "posterior_sd = np.sqrt(posterior_var)\n",
    "plt.axvline(posterior_mean - posterior_sd, linestyle='--', color='g', label='posterior mean - 1 s.d.')\n",
    "plt.axvline(posterior_mean + posterior_sd, linestyle='--', color='g', label='posterior mean + 1 s.d.')\n",
    "plt.legend(loc='upper left');\n",
    "print('posterior mean:', posterior_mean)\n",
    "print('posterior s.d.:', posterior_sd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the results of Bayesian estimation are sensitive, to varying degree in each specific case, to the choice of prior distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After all, according to Stephen Senn's *Statistical basis of public policy &mdash; present remembrance of priors past is not the same as a true prior*, British Medical Journal, 1997, **\"*A Bayesian is one who, vaguely expecting a horse, and catching a glimpse of a donkey, strongly concludes he has seen a mule.*\"**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sequential Bayesian updates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous section we saw that, starting with the prior\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha, \\beta\\right),$$\n",
    "we end up with another Beta-distributed posterior,\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha + \\sum x_i, \\beta + n - \\sum x_i\\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What would happen if, instead of observing all ten coin tosses at once, we considered each coin toss in turn, obtained our posterior, and used that posterior as a prior for an update based on the information from the next coin toss?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above two formulae give the answer to this question. We start with our initial prior,\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha, \\beta\\right),$$\n",
    "then, substituting $n = 1$ into the second formula, we get\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha + x_1, \\beta + 1 - x_1\\right).$$\n",
    "Using this posterior as a prior before the second coin toss, we obtain the next posterior as\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha + x_1 + x_2, \\beta + 2 - x_1 - x_2\\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Proceeding along these lines, after all ten coin tosses, we end up with\n",
    "$$\\text{Beta}\\left(\\theta \\, | \\, \\alpha + \\sum x_i, \\beta + n - \\sum x_i\\right),$$\n",
    "the same result that we would got if we processed all ten coin tosses as a single \"batch\" as we did the previous section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This insight forms the basis for a **sequential** or **iterative** application of Bayes's theorem, sequential Bayesian updates, the foundation of real-time **Bayesian filtering**."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
